{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scrublet as scr\n",
    "import scipy.io\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Download 8k PBMC data set from 10X Genomics\n",
    "(Only need to do this once)\n",
    "\n",
    "Download raw data by navigating to the following URL in your web browser:\n",
    "http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_filtered_gene_bc_matrices.tar.gz\n",
    "\n",
    "Or use wget:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-07-16 16:49:19--  http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_filtered_gene_bc_matrices.tar.gz\n",
      "Resolving cf.10xgenomics.com (cf.10xgenomics.com)... 13.33.35.84, 13.33.35.97, 13.33.35.175, ...\n",
      "Connecting to cf.10xgenomics.com (cf.10xgenomics.com)|13.33.35.84|:80... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 37558165 (36M) [application/x-tar]\n",
      "Saving to: ‘pbmc8k_filtered_gene_bc_matrices.tar.gz’\n",
      "\n",
      "pbmc8k_filtered_gen 100%[===================>]  35.82M  7.21MB/s    in 6.4s    \n",
      "\n",
      "2018-07-16 16:49:25 (5.60 MB/s) - ‘pbmc8k_filtered_gene_bc_matrices.tar.gz’ saved [37558165/37558165]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_filtered_gene_bc_matrices.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncompress:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar xfz pbmc8k_filtered_gene_bc_matrices.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load counts matrix and gene list\n",
    "The first time this is run, the counts matrix is loaded from the mtx file. An npz file is saved for fast loading in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expression matrix shape: 8381 rows, 33694 columns\n",
      "Number of genes in gene list: 33694\n"
     ]
    }
   ],
   "source": [
    "input_dir = 'filtered_gene_bc_matrices/GRCh38/'\n",
    "\n",
    "# The raw counts matrix (E) should be a scipy sparse CSC matrix\n",
    "# with cells as rows and genes as columns\n",
    "\n",
    "if os.path.isfile(input_dir + '/matrix.npz'):\n",
    "    E = scipy.sparse.load_npz(input_dir + '/matrix.npz')\n",
    "else:\n",
    "    E = scipy.io.mmread(input_dir + '/matrix.mtx').T.tocsc()\n",
    "    scipy.sparse.save_npz(input_dir + '/matrix.npz', E, compressed=True)\n",
    "\n",
    "genes = np.array(scr.load_genes(input_dir + 'genes.tsv', delimiter='\\t', column=1))\n",
    "\n",
    "print('Expression matrix shape: {} rows, {} columns'.format(E.shape[0], E.shape[1]))\n",
    "print('Number of genes in gene list: {}'.format(len(genes)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check that the distribution of total counts per cell looks reasonable (i.e., background has been filtered out)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Number of cells')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAADXCAYAAACu9hJ0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFPZJREFUeJzt3X2UXVV5x/Hvj4Q3FQghgUXz0kSJL9RXnAJqVTQoBJAghRZUjJhltIJSsZVgralFa/ANoSqaQiQoC8RIJUgUYiSAXQaSIIaXoEwxJUOQJCZAAAEDT/84e+AS7tw5c+fu+zLz+6w1656zzz7nPMlZ88y+++yzjyICMzNrrB1aHYCZ2VDk5GpmloGTq5lZBk6uZmYZOLmamWXg5GpmloGTq5lZBk6uZmYZOLmamWUwstUB5DBmzJiYNGlSq8MwsyFm1apVmyJibJm6QzK5Tpo0iZUrV7Y6DDMbYiT9X9m67hYwM8sgW3KVNF/SBkm3V5R9WdJdklZL+m9Joyq2nSmpW9JvJR1WUX54KuuWNDtXvGZmjZSz5XoRcPh2ZUuAV0bEq4HfAWcCSNofOAH4q7TPtySNkDQC+CYwDdgfODHVNTNra9mSa0TcAGzeruzaiNiWVpcD49PydOCyiHgiIn4PdAMHpp/uiLgnIp4ELkt1zczaWiv7XD8I/DQtjwPWVWzrSWV9lZuZtbWWJFdJ/wJsAy7pLapSLWqUVzvmLEkrJa3cuHFjYwI1M6tT04diSZoBHAVMjWdfg9ADTKioNh5Yn5b7Kn+OiJgHzAPo6ury6xUabNLsq6uWr517ZJMjMesMTW25SjocOAM4OiIeq9i0CDhB0s6SJgNTgJuBFcAUSZMl7URx02tRM2M2M6tHtparpEuBQ4AxknqAORSjA3YGlkgCWB4RH4mIOyRdDtxJ0V1wSkQ8lY5zKnANMAKYHxF35IrZzKxRsiXXiDixSvGFNep/AfhClfLFwOIGhmZmlp2f0DIzy2BIzi1gzdPXjS7wzS4b3txyNTPLwMnVzCwDJ1czswycXM3MMnByNTPLwMnVzCwDJ1czswycXM3MMnByNTPLwMnVzCwDJ1czswycXM3MMnByNTPLwLNi2TNqzXBlZgPjlquZWQbZkquk+ZI2SLq9omy0pCWS7k6fe6ZySTpPUrek1ZIOqNhnRqp/d3q5oZlZ28vZcr0IOHy7stnA0oiYAixN6wDTKF5KOAWYBZwPRTKmePfWQcCBwJzehGxm1s6yJdeIuAHYvF3xdGBBWl4AHFNRfnEUlgOjJO0LHAYsiYjNEbEFWMLzE7aZWdtpdp/rPhFxP0D63DuVjwPWVdTrSWV9lZuZtbV2uaGlKmVRo/z5B5BmSVopaeXGjRsbGpyZ2UA1O7k+kL7ukz43pPIeYEJFvfHA+hrlzxMR8yKiKyK6xo4d2/DAzcwGotnJdRHQe8d/BnBlRfn706iBg4GHUrfBNcA7Je2ZbmS9M5WZmbW1bA8RSLoUOAQYI6mH4q7/XOBySTOBe4HjU/XFwBFAN/AYcDJARGyWdBawItX794jY/iaZmVnbyZZcI+LEPjZNrVI3gFP6OM58YH4DQzMzy65dbmiZmQ0pTq5mZhk4uZqZZeDkamaWgZOrmVkGTq5mZhn0m1wlfUnS7pJ2lLRU0iZJ72tGcGZmnapMy/WdEfEwcBTF46gvBf45a1RmZh2uTHLdMX0eAVzqJ6TMzPpX5gmtqyTdBfwJ+KikscDjecMyM+ts/bZcI2I28AagKyL+DDxKMbm1mZn1oc+Wq6Rjq5RVrl6RIyAzs6GgVrfAu2psC5xczcz61GdyjYiTmxmIDT2TZl9dtXzt3CObHIlZ89XqFji91o4R8bXGh2NmNjTU6hbYrWlRmJkNMbW6BT7XzEDMzIaSMo+/vjQ99np7Wn+1pM8M5qSSPiHpDkm3S7pU0i6SJku6SdLdkn4gaadUd+e03p22TxrMuc3MmqHME1r/BZwJ/BkgIlYDJ9R7QknjgI9TjJt9JTAiHe9s4JyImAJsAWamXWYCWyJiP+CcVM/MrK2VSa4viIibtyvbNsjzjgR2lTQSeAFwP/B2YGHavgA4Ji1PT+uk7VO13YBbM7N2U+bx102SXkIxthVJx1Ekw7pExH2SvkLx9tc/AdcCq4AHI6I3afcA49LyOGBd2nebpIeAvYBN9cYw3PU1RMrMGqdMcj0FmAe8XNJ9wO+BuqcclLQnRWt0MvAg8ENgWpWq0btLjW2Vx50FzAKYOHFiveGZmTVEv8k1Iu4BDpX0QmCHiNg6yHMeCvw+IjYCSLoCeCMwStLI1HodD6xP9XuACUBP6kbYA3jezFwRMY/ijwBdXV3PS75mZs1UZrTAf0gaFRGPRsRWSXtK+vwgznkvcLCkF6S+06nAncB1wHGpzgzgyrS8KK2Ttv8iIpw8zaytlbmhNS0iHuxdiYgtFHO71iUibqK4MXULcFuKYR5wBnC6pG6KPtUL0y4XAnul8tOB2fWe28ysWcr0uY6QtHNEPAEgaVdg58GcNCLmAHO2K74HOLBK3ceB4wdzPjOzZiuTXL8PLJX0XYobSR/k2aFRZmZWRZkbWl+StJriRpSAsyLimuyRmZl1sDItVyLiZ8DPMsdiZjZklLmhZWZmA+TkamaWQZ/JVdLS9OmJUszMBqhWn+u+kt4KHC3pMrZ7DDUibskamZlZB6uVXD9LMWB/PLD9K12CYhYrMzOrotabCBYCCyX9a0Sc1cSYzMw6XplxrmdJOhp4SypaFhE/yRuWmVlnKzNxyxeB0ygmV7kTOC2VmZlZH8o8RHAk8NqIeBpA0gLg1xSvfjEzsyrKjnMdVbG8R45AzMyGkjIt1y8Cv5Z0HcVwrLfgVquZWU1lbmhdKmkZ8NcUyfWMiPhD7sDMzDpZ2Ylb7qd4I4CZmZXguQXMzDJwcjUzy6BmcpW0g6TbG31SSaMkLZR0l6Q1kt4gabSkJZLuTp97prqSdJ6kbkmrJR3Q6HjMzBqtZnJNY1t/I2lig897LvCziHg58BpgDcU8BksjYgqwlGdfRDgNmJJ+ZgHnNzgWM7OGK3NDa1/gDkk3A4/2FkbE0fWcUNLuFMO5PpCO8yTwpKTpwCGp2gJgGcUbYacDF6fXaS9Prd590002M7O2VCa5fq7B53wxsBH4rqTXAKsoHq/dpzdhRsT9kvZO9ccB6yr270llz0mukmZRtGyZOLHRDW0zs4Hp94ZWRFwPrAV2TMsrgMHM5ToSOAA4PyJeR9Eanl2jvqqURZU450VEV0R0jR07dhDhmZkNXr8tV0kfomgRjgZeQtFq/DYwtc5z9gA9EXFTWl9IkVwf6P26L2lfYENF/QkV+48H1td5bmsDk2Zf3ee2tXOPbGIkZvmUGYp1CvAm4GGAiLgb2LvmHjWkp7vWSXpZKppKMdvWImBGKpsBXJmWFwHvT6MGDgYecn+rmbW7Mn2uT0TEk1Lx7VzSSKp8LR+gjwGXSNoJuAc4mSLRXy5pJnAvcHyquxg4AugGHkt1zczaWpnker2kTwO7SnoH8FHgqsGcNCJuBbqqbHpeV0MaJXDKYM5nZtZsZboFZlPc3b8N+DBFS/IzOYMyM+t0ZWbFejpNkH0TRXfAb1Nr0szM+lBmtMCRFKMD/pdiWNRkSR+OiJ/mDs7MrFOV6XP9KvC2iOgGkPQS4GrAydXMrA9l+lw39CbW5B6eHYNqZmZV9NlylXRsWrxD0mLgcoo+1+MpntKyNlZroL6Z5VerW+BdFcsPAG9NyxuBPbNFZGY2BPSZXCPCg/XNzOpUZrTAZIonqiZV1q93ykEzs+GgzGiBHwMXUjyV9XTecMzMhoYyyfXxiDgveyRmZkNImeR6rqQ5wLXAE72FETGYOV3NzIa0Msn1VcBJwNt5tlsg0rqZmVVRJrm+G3hxeteVmZmVUOYJrd8Ao3IHYmY2lJRpue4D3CVpBc/tc/VQLDOzPpRJrnOyR2FmNsSUmc/1+hwnljQCWAncFxFHpYcVLqN4EeItwEnp9TI7AxcDrwf+CPx9RKzNEZOZWaP02+cqaaukh9PP45KekvRwA859GrCmYv1s4JyImAJsAWam8pnAlojYDzgn1TMza2v9JteI2C0idk8/uwB/C3xjMCeVNB44ErggrYtiaNfCVGUBcExanp7WSdunqvdtiWZmbapMn+tzRMSPJc0e5Hm/DnwK2C2t7wU8GBHb0noPMC4tjwPWpXNvk/RQqr+p8oCSZgGzACZOnDjI8KxV+poqce3cI5scidnglJm45diK1R0o3tpa9zu0JB1FMQH3KkmH9BZXqRoltj1bEDEPmAfQ1dXld3yZWUuVablWzuu6DVhL8VW9Xm8CjpZ0BLALsDtFS3aUpJGp9ToeWJ/q9wATgB5JI4E9gM2DOL+ZWXZlRgs0dF7XiDgTOBMgtVz/KSLeK+mHwHEUIwZmAFemXRal9V+l7b/w22fNrN3Ves3LZ2vsFxFxVoNjOQO4TNLngV9TTHNI+vyepG6KFusJDT6vmVnD1Wq5Plql7IUUQ6P2AgadXCNiGbAsLd8DHFilzuMU7+0yM+sYtV7z8tXeZUm7UYxLPZnia/tX+9rPzMz66XOVNBo4HXgvxVjTAyJiSzMCMzPrZLX6XL8MHEsxvOlVEfFI06IyM+twtZ7Q+iTwF8BngPUVj8BubdDjr2ZmQ1atPtcyc71ai/X1RJOZtdaAH381a4Vaf0T8aKy1I7dOzcwycHI1M8vAydXMLAMnVzOzDJxczcwycHI1M8vAydXMLAMnVzOzDJxczcwycHI1M8vAj79ax/MbY60dNb3lKmmCpOskrZF0h6TTUvloSUsk3Z0+90zlknSepG5JqyUd0OyYzcwGqhXdAtuAT0bEK4CDgVMk7Q/MBpZGxBRgaVoHmAZMST+zgPObH7KZ2cA0PblGxP0RcUta3gqsAcZRvK57Qaq2ADgmLU8HLo7CcopXcO/b5LDNzAakpTe0JE0CXgfcBOwTEfdDkYCBvVO1ccC6it16Utn2x5olaaWklRs3bswZtplZv1qWXCW9CPgR8I8RUevNBqpSFs8riJgXEV0R0TV27NhGhWlmVpeWJFdJO1Ik1ksi4opU/EDv1/30uSGV9wATKnYfD6xvVqxmZvVo+lAsSQIuBNZExNcqNi0CZgBz0+eVFeWnSroMOAh4qLf7wKwWv73AWqkV41zfBJwE3Cbp1lT2aYqkermkmcC9wPFp22LgCKAbeAw4ubnhmpkNXNOTa0T8kur9qABTq9QP4JSsQZmZNZgffzUzy8DJ1cwsA88t0AFq3Zgxs/bklquZWQZOrmZmGTi5mpll4D7XNuK+VbOhw8m1yZxA24Mn2Lbc3C1gZpaBW65mFTwfgTWKk2sm/vo/9LgrwQbC3QJmZhk4uZqZZeBuAbNBcj+tVeOWq5lZBk6uZmYZOLmamWXQMX2ukg4HzgVGABdExNwWh2TWLw/fGr46IrlKGgF8E3gHxdtgV0haFBF3tjYys/r4JtjQ1xHJFTgQ6I6IewDSm2CnAy1Nrn5QwHJwa3do6JTkOg5YV7HeQ/Ga7WdImgXMSquPSPoD8FCVY+3RR/kYYNPgQ224vuJt9bHr2bfsPv3Vq3d7R197nZ3nuE3avxnXvp5tA732f1m6ZkS0/Q/Fa7YvqFg/CfjPfvaZN8Dyla3+dw4k3lYfu559y+7TX716t/vaN+a47Xrt69mW89p3ymiBHmBCxfp4YH0/+1w1wPJ2lTPewRy7nn3L7tNfvXq3+9o35rjteu3r3ZaFUvZua5JGAr8DpgL3ASuA90TEHQ08x8qI6GrU8axz+NoPXzmvfUf0uUbENkmnAtdQDMWa38jEmsxr8PGsc/jaD1/Zrn1HtFzNzDpNp/S5mpl1FCdXM7MMnFzNzDJwcjUzy8DJtQ+SXiHp25IWSvqHVsdjzSXphZJWSTqq1bFYc0g6RNKN6ff+kMEeb1glV0nzJW2QdPt25YdL+q2kbkmzASJiTUR8BPg7wGMgO9xArn1yBnB5c6O0RhvgdQ/gEWAXigeXBmVYJVfgIuDwyoKKGbemAfsDJ0raP207GvglsLS5YVoGF1Hy2ks6lGJSoAeaHaQ13EWU/52/MSKmUfxh/dxgTzyskmtE3ABs3q74mRm3IuJJoHfGLSJiUUS8EXhvcyO1RhvgtX8bcDDwHuBDkobV78lQMpDrHhFPp+1bgJ0He+6OeEIrs6ozbqU+l2Mp/pMXtyAuy6/qtY+IUwEkfQDYVPFLZ0NDX7/zxwKHAaOAbwz2JE6uoCplERHLgGXNDcWarOq1f2Yh4qLmhWJN1Nfv/BXAFY06ib/u1Dfjlg0NvvbDU1Ouu5NrMcPWFEmTJe0EnAAsanFM1hy+9sNTU677sEquki4FfgW8TFKPpJkRsQ3onXFrDXB5hhm3rMV87YenVl53z4plZpbBsGq5mpk1i5OrmVkGTq5mZhk4uZqZZeDkamaWgZOrmVkGTq7WFiTtJenW9PMHSfdVrO9Upf5oSR8pcdyRkh7ME3W/5z5d0i6tOLe1nse5WtuR9G/AIxHxlRp19gMWRsRr+znWSIrJV0Y1Nsr+SeoBXhkRLUnu1lpuuVrbk/QpSbenn4+l4rkUT93cKmmupN0l/ULSLZJWl3mDgKSTU93fSPpuKpss6bpUvkTS+FT+fUnHVOz7SPo8VNJSSVekyZcvTuWfAPYGbpT089SC/p6k29K/4+ON/V+yduNZsaytSTqQYj7dA4ERwM2SrgdmA/v1tlwl7UgxJ+dWSXsD/wP8pMZxX0MxKfIbI2KzpNFp07eACyLiEkmzgK8Dx/UT5gEUky5vAJZLOjgizpH0SeDNEfGgpIOAMRHxqnT+prekrbnccrV292bgRxHxWERsBX4M/E2VegLOlrQauBaYIGlMjeO+HfhBRGwG6P0EDqKYPBng4nT+/iyPiPsj4ingVmBSlTrdFC3tcyUdBjxU4rjWwZxcrd1Vm3uzmvcDewAHpNbsJop3IdU67kBuOGwj/b6k14RUfut7omL5Kap8I4yIPwKvpnht0MeB7wzg3NaBnFyt3d0AvFvSrpJeRPEalhuBrcBuFfX2ADZExDZJ76CYbb6WnwMn9HYHVHQLLKd4KSXA+9L5AdYCr0/L76booujPMzFKGktxA/mHwByKrgQbwtznam0tIm5O08atSEXnR8RtAJJWSroNuBr4GnCVpJXALcDd/Rx3taQvATdI2gasAmZSTEV3oaQzKV5QeHLa5TvAlSlxX8tzW6t9mQf8XNI64FPpuL0t5jPK/Q9Yp/JQLDOzDNwtYGaWgZOrmVkGTq5mZhk4uZqZZeDkamaWgZOrmVkGTq5mZhn8P03I/hKRUzg6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "total_counts = E.sum(1).A.squeeze()\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (5, 3))\n",
    "ax.hist(total_counts, bins = np.logspace(3, 5, 40))\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('Total counts')\n",
    "ax.set_ylabel('Number of cells')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate doublet scores\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulating doublets\n",
      "Total counts normalizing\n",
      "Finding highly variable genes\n",
      "Filtering genes from 33694 to 1697\n",
      "Applying z-score normalization\n",
      "Running PCA\n",
      "Building kNN graph and calculating doublet scores\n",
      "Elapsed time: 19.0 seconds\n"
     ]
    }
   ],
   "source": [
    "# filtering/preprocessing parameters:\n",
    "min_counts = 2\n",
    "min_cells = 3\n",
    "vscore_percentile = 85\n",
    "n_pc = 30\n",
    "\n",
    "# doublet detector parameters:\n",
    "expected_doublet_rate = 0.06 \n",
    "sim_doublet_ratio = 3\n",
    "n_neighbors = 50\n",
    "\n",
    "t0 = time.time()\n",
    "\n",
    "scrublet_results = scr.compute_doublet_scores(\n",
    "    E, \n",
    "    min_counts = min_counts, \n",
    "    min_cells = min_cells, \n",
    "    vscore_percentile = vscore_percentile, \n",
    "    n_prin_comps = n_pc,\n",
    "    scaling_method = 'zscore',\n",
    "    expected_doublet_rate = expected_doublet_rate,\n",
    "    sim_doublet_ratio = sim_doublet_ratio,\n",
    "    n_neighbors = n_neighbors, \n",
    "    use_approx_neighbors = True, \n",
    "    get_doublet_neighbor_parents = False\n",
    ")\n",
    "\n",
    "\n",
    "t1 = time.time()\n",
    "print('Elapsed time: {:.1f} seconds'.format(t1 - t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get UMAP embedding to help visualize the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# UMAP: https://github.com/lmcinnes/umap\n",
    "import umap\n",
    "\n",
    "embedding = umap.UMAP(n_neighbors=10).fit_transform(scrublet_results['pca_observed_cells'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set doublet score threshold and visualize results\n",
    "To call doublets, manually set a threshold between the two peaks of the simulated doublet histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "349/8381 = 4.2% of cells are predicted doublets.\n",
      "60.0% of doublets are predicted to be detectable.\n",
      "Predicted overall doublet rate = 6.9%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmYHWWZ///3p093J52tE7KwJSFgIBK2AA2EURFnUBZB5lIcQWRkBkR0UByXGRxmZFF/o/MdcUNlgjBBZR10FBBEVBBxIJBggAQEAgQSErKSTtJJer1/f1R1ctLp5aT7bH3687quc+Wcqqeq7nOSunPXU09VKSIwMzMzqyRVpQ7AzMzMLN9c4JiZmVnFcYFjZmZmFccFjpmZmVUcFzhmZmZWcVzgmJmZWcVxgWPdknSlpJ+UOo7dIel8SY/keZ0nSlqe9XmppJPyuQ2zciXpXEm/LtC650r6SiHW3c22dis35Gs/7+s7SgpJ0we6HeueC5whKt3hn5G0RdIbkn4gaWyp4zKz4pL0dkn/J6lR0npJf5R0DEBE3BwR7ymDGB+SdGGp4yiVrgdalhsXOEOQpM8BXwe+ANQDs4H9gAck1RYxjupibcvMdiVpDHAP8F1gD2Bf4CqguZRxmeWDC5whJk1oVwGfiohfRURrRCwF/oakyPlIVvPhkm6XtEnSk5KOyFrPP0t6PZ33vKS/SqdXSbpM0kuS1km6Q9Ie6bxpaZfsBZJeA34n6VeSLukS41OS3p++f6ukB9Ijy+cl/U1Wu/GS7pK0UdLjwFv6+O6dR6obJC2TdH46fZik/5T0mqRVkq6TVJfDb3mspPnp9ldJuqavZczKzEEAEXFrRLRHxNaI+HVEPA27ntpJ999PSnox3fe/LOktkh5N94M7Og+Sujst1NMpGUnjJN0jaY2kN9P3k9N5XwXeAVwrabOka9Pp+cwN50l6Nc1Zl3eZN0zStyStSF/fkjRsN77jhDTOTZJ+L2m/HmLoNg9JGgncB+yTfv/NkvZx/umbC5yh5y+A4cDPsidGxGaSnejdWZPPBP6H5MjuFuDnkmokzQAuAY6JiNHAycDSdJlPA38NvBPYB3gT+F6XGN4JHJwudwtwTucMSTNJCq1fpjv2A2mbSWm770s6JG3+PWAbsDfw9+mrW5Kmpt/vu8BEYBawMJ39dZJEPwuYTnIU+6We1pXl28C3I2IMSQK9I4dlzMrJC0C7pJsknSppXA7LnAIcTdLz+0/AHOBcYApwKFn7826oAv6bZN+fCmwFrgWIiMuBPwCXRMSoiLgkz7lhJvAD4DySnDUemJzV5PL0u84CjgCOBf51N77bucCXgQkkOefmHtp1m4ciogk4FViRfv9REbEC558+ucAZeiYAayOirZt5K9P5nRZExJ0R0QpcQ1IYzQbagWHATEk1EbE0Il5Kl/k4cHlELI+IZuBK4CztfDrqyohoioitwP8Cs7KOas4FfpYuezqwNCL+OyLaIuJJ4Kfp+jLAB0gTQEQsAm7q5XufC/wmPVJtjYh1EbFQkoCPAf8YEesjYhPw/wFn9/lLQiswXdKEiNgcEY/lsIxZ2YiIjcDbgQCuB9akPR979rLY1yNiY0QsBhYBv46IlyOikeQg4sh+xLEuIn4aEVvSffCrJAdCPclnbjgLuCciHk7zzr8BHVnzzwWujojVEbGGpAf8vN34er/MWvflwPGSpmQ36Gcecv7pgwucoWctSZdpd+Nf9k7nd1rW+SYiOoDlwD4RsQT4DEnxslrSbZL2SZvuB/xvehpoA/AcSUG0Zw/r3QT8kh078tnsOMLZDziuc13p+s4F9iLphanOXhfwai/fewrwUjfTJwIjgAVZ2/hVOr0vF5Accf1Z0hOSTs9hGbOyEhHPRcT5ETGZpAdmH+BbvSyyKuv91m4+j9rdGCSNkPRf6WmijcDDwNi0WOlOPnPDPuyck5qAdV3mZy//ajotV9nr3gys72b5/uQh558+uMAZeh4lGUD4/uyJaZfvqcBvsyZPyZpfRdJtuwIgIm6JiLeTJJog6V6FZGc+NSLGZr2GR8TrWevt+gj7W4FzJB0P1AEPZq3r913WNSoiPgGsAdqyYyTp2u7JMro/D7+WJCkfkrWN+ojoM0lHxIsRcQ5JF/nXgTvT39FsUIqIPwNzSQqdgWoi+U8bAEl79dL2c8AM4Lj0lMsJnYt1htalfT5zw0p2znUjSE5TdVpBkuey17UifZ/Ld8xe9yiSU/4rurTpKw91/f7OPzlwgTPEpN3IVwHflXRKOqZmGslYm+XAj7OaHy3p/Wlvz2dICqPHJM2Q9JfpQLttJDtme7rMdcBXO085SZoo6cw+wrqXJIFcDdye9hZBcnXHQekAwJr0dYykgyOinWQc0ZXp0d9M4KO9bONm4CRJfyOpOh2EOCvd1vXANyVNSmPeV9LJfcSMpI9ImpiuY0M6ub23ZczKiZKBup/TjgG9U0jGs+TjdMdTwCGSZkkaTtLj25PRJHlkg5KLEq7oMn8VcEDW53zmhjuB05VchFBLkoey/2+8FfjXNJdNIBmf13mPsFy+42lZ6/4yMC8isnuXyCEPrQLGS6rvXMb5p28ucIagiPgP4F+A/wQ2AvNIjoj+Kj1P3OkXwIdIBgqfB7w/HY8zDPgayVHHGyRHEP+SLvNt4C7g15I2kSTK4/qIp5kkIZ1EMmiwc/om4D0kp61WpNv6erp9SAY6j0qnzyUZpNjTNl4DTiM5UlxPMtiv86qwfwaWkBRvG4HfkBxN9uUUYLGkzen3PjsituWwnFm52ESyf86T1ESyvy4i2U8GJCJeICkWfgO8CPR2o71vkfTerk1j+FWX+d8mGV/zpqTv5Dk3LAb+gST3rCTJd9n3nPkKMB94GngGeDKdlut3vIWkYFtPMjj73B5C6TEPpT1rtwIvp6ew9sH5p0+K2KXny8zMzGxQcw+OmZmZVRwXOGZmZlZxXOCYmZlZxXGBY2ZmZhWnIh92OGHChJg2bVqpwyiJBQsWcPTRR5c6DLMeLViwYG1E5HIjxZJxDnEOsfKVaw6pyKuoGhoaYv78+aUOoyQkUYl/p1Y5JC2IiIZSx9Eb5xDnECtfueYQn6IyMzOziuMCx8zMzCqOCxwzMzOrOBVV4Eg6Q9KcxsbGUodiZmZmJVRRBU5E3B0RF9XX1/fd2MzMzCpWRRU4ZmZmZlCh98EZiKuuumqnz1dccUWJIjEzs4FwPh/a3INjZmZmFccFjpmZmVUcFzhmZmZWcVzgmJmZWcUZ0oOMuw5AMzMzs8rgHhwzMzOrOEO6B8fMzCybLy2vHO7BMbNBS9JwSY9LekrSYkm7nHeWdL6kNZIWpq8LSxGrmRWXe3DMbDBrBv4yIjZLqgEekXRfRDzWpd3tEXFJCeIzsxIp+x4cSQdIukHSnaWOxczKSyQ2px9r0leUMCQzKxMlKXAk3ShptaRFXaafIul5SUskXQYQES9HxAWliNPMyp+kjKSFwGrggYiY102zD0h6WtKdkqYUOUQzK4FS9eDMBU7JniApA3wPOBWYCZwjaWbxQzOzwSQi2iNiFjAZOFbSoV2a3A1Mi4jDgd8AN3W3HkkXSZovaf6aNWsKG7SZFVxJCpyIeBhY32XyscCStMemBbgNODPXdTo5mQ1tEbEBeIguB08RsS4imtOP1wNH97D8nIhoiIiGiRMnFjRWMyu8chpkvC+wLOvzcuA4SeOBrwJHSvpiRPx7dwtHxBxgDkBDQ4PPwZsNAZImAq0RsUFSHXAS8PUubfaOiJXpx/cBzxU5TMsDX75tu6ucChx1My0iYh1wcbGDMbNBYW/gpvQUdxVwR0TcI+lqYH5E3AV8WtL7gDaSnuPzSxatmRVNORU4y4HswX+TgRW7swJJZwBnTJ8+PZ9xmVmZioingSO7mf6lrPdfBL5YzLjMrPTK6TLxJ4ADJe0vqRY4G7hrd1YQEXdHxEX19fUFCdDMzMwGh1JdJn4r8CgwQ9JySRdERBtwCXA/yTnyOyJi8W6u9wxJcxobG/MftJmZmQ0aJTlFFRHn9DD9XuDeAaz3buDuhoaGj/V3HV15YJuZ7Q7nDLPyUE6nqMzMzMzywgWOmZmZVZxyuopqwHwVlZmZlZuupy3Bpy6LoaJ6cHwVlZmZmUGFFThmZmZmUGEFji8TNzMzM6iwAsenqMzMzAwqrMAxMzMzAxc4ZmZmVoFc4JjZoCVpuKTHJT0labGkXa7HlTRM0u2SlkiaJ2la8SM1s2KrqALHg4zNhpxm4C8j4ghgFnCKpNld2lwAvBkR04FvAl8vcoxmVgIVVeB4kLHZ0BKJzenHmvQVXZqdCdyUvr8T+CtJKlKIZlYiFVXgmNnQIykjaSGwGnggIuZ1abIvsAwgItqARmB8caM0s2KrqEc1mNnQExHtwCxJY4H/lXRoRCzKatJdb03XXh4kXQRcBDB16tSCxGo981PYLd/cg2NmFSEiNgAPAad0mbUcmAIgqRqoB9Z3s/yciGiIiIaJEycWOFozK7SKKnA8yNhsaJE0Me25QVIdcBLw5y7N7gI+mr4/C/hdROzSg2NmlaWiChwPMjYbcvYGHpT0NPAEyRiceyRdLel9aZsbgPGSlgCfBS4rUaxmVkQeg2Nmg1ZEPA0c2c30L2W93wZ8sJhxmVnpucDZTR4IZ2ZmVv4q6hSVmZmZGbgHx8ysoLr2+oJ7fs2KwT04ZmZmVnEqqgdH0hnAGdOnTy91KGZmZiU3lMeNVlQPji8TNzMzM8ihwJG0RzECMbOhy3nGzPItlx6ceZL+R9JpfgKvmRWI84yZ5VUuBc5BwBzgPGCJpP9P0kGFDcvMhhjnGTPLqz4LnEg8EBHnABeSPNPlcUm/l3R8wSM0s4rnPGNm+dbnVVSSxgMfITmyWgV8iuThdbOA/wH2L2SAZlb5nGfMLN9yuUz8UeDHwF9HxPKs6fMlXVeYsMxsiHGeMbO8ymUMzr9GxJezk46kDwJExNcLFpmZDSX9yjOSpkh6UNJzkhZLurSbNidKapS0MH19qbt1mVllyaXAuaybaV/MdyD5IOkMSXMaGxtLHYqZ7Z7+5pk24HMRcTAwG/gHSTO7afeHiJiVvq4eSKBmNjj0eIpK0qnAacC+kr6TNWsMSVIpOxFxN3B3Q0PDx0odi5n1baB5JiJWAivT95skPQfsCzxbgHDNbBDpbQzOCmA+8D5gQdb0TcA/FjIoMxsy8pZnJE0DjgTmdTP7eElPpdv7fEQs7k+wZjZ49FjgRMRTwFOSbo6IsuyxMbPBLV95RtIo4KfAZyJiY5fZTwL7RcRmSacBPwcO7GYdFwEXAUydOrW/oZhZmejtFNUdEfE3wJ8kRfYskttWHF7w6AaBrg8yg6H1MDOzgchHnpFUQ1Lc3BwRP+s6P7vgiYh7JX1f0oSIWNul3RySmw3S0NAQmNmg1tspqs6rEU4vRiBmNiQNKM+kj3W4AXguIq7poc1ewKqICEnHklxcsa4/2zOzwaO3U1Qr07drga0R0ZHeOv2twH3FCM7MKlse8szbSG4O+Iykhem0fwGmpuu/DjgL+ISkNmArcHZEuIfGrMLlcqO/h4F3SBoH/JZkQOCHgHMLGZiZDSn9yjMR8QjJ6aze2lwLXJunOM1skMilwFFEbJF0AfDdiPgPSX8qdGBmNqQMqTzTdexepY/bG2rf18pDLjf6U/qwu3OBX6bTcimMzMxy5TxjZnmVS4FzKckdRf83IhZLOgB4sLBhmdkQ4zxjZnnV5xFSRDxMcn688/PLwKcLGZSZDS3OM2aWb30WOOkVDZ8HpmW3j4i/LFxYZjaUOM+YWb7lco77f4DrgB8C7YUNZ1eSRgLfB1qAhyLi5mLHYGYFV9I8Y2aVJ5cCpy0ifpDPjUq6keTGXqsj4tCs6acA3wYywA8j4mvA+4E7I+JuSbcDLnDMKk/e84yZDW25DDK+W9InJe0taY/O1wC3Oxc4JXuCpAzwPeBUYCZwjqSZwGRgWdrMR3ZmlakQecbMhrBcenA+mv75haxpARzQ341GxMPpk3+zHQssSQcXIuk24ExgOUmRs5DcCjIzG3zynmfMrHwV4zmOuVxFtX9et9izfdnRUwNJYXMc8B3gWknvBe7uaWE/Cdhs8CpinilLvhGeWf7lchXVCOCzwNSIuEjSgcCMiLgnz7F0d7v1iIgm4O/6WthPAjYbvIqYZ8xsiMjllM9/k1zB9Bfp5+XAVwoQy3JgStbnycCK3VmBpDMkzWlsbMxrYGZWcMXKM2Y2RORS4LwlIv4DaAWIiK308XC7fnoCOFDS/pJqgbOBu3ZnBRFxd0RcVF9fX4DwzKyAipVnzGyIyGWQcYukOpIBf0h6C9A8kI1KuhU4EZggaTlwRUTcIOkS4H6Sy8RvjIjFA9lOqfh8utlu61eekTQF+BGwF9ABzImIb3dpI5LbT5wGbAHOj4gn8xu+mZWbXAqcK4BfAVMk3Qy8DTh/IBuNiHN6mH4vcG9/1yvpDOCM6dOn93cVZlYa/c0zbcDnIuJJSaOBBZIeiIhns9qcChyYvo4DfpD+aWYVLJerqB6Q9CQwm6TL+NKIWFvwyPohIu4G7m5oaPhYqWMxs9z1N89ExEpgZfp+k6TnSK7IzC5wzgR+FBEBPCZprKS902XNrEL1WOBIOqrLpM5kMFXSVHfxmtlA5TPPpPfWOhKY12VWd7eg2DdrW53L+1YTZhWktx6cb6R/DgcagKdIjqwOJ0kgby9saLvPp6jMBp285BlJo4CfAp+JiI1dZ3ezyC63kvCtJswqS49XUUXEuyLiXcCrwFER0RARR5McIS0pVoC7w1dRmQ0u+cgzkmpIipubI+Jn3TQZ8C0ozGzwyeUy8bdGxDOdHyJiETCrcCGZ2RDUrzyTXiF1A/BcRFzTQ7O7gL9VYjbQ6PE3ZpUvl6uonpP0Q+AnJN26HwGeK2hUZjbU9DfPvA04D3hG0sJ02r8AUwEi4jqSKzNPI+kR2kIOd0Y3s8EvlwLn74BPAJemnx8mucyy7HgMjtmg1a88ExGP0McNAdOrp/5hoAGa2eCSy2Xi24Bvpq+y5svEzQanwZRnbGfFeCq0WX/k0oNjA+Q7G5uZmRVXLoOMzczMzAaViipw/DRxMzMzg34WOOkdP8uO74NjVjnKNc+Y2eDQ3zE4vV61YL3zoDyznDjPmFm/9asHJyL+K9+BmJllc54xs4Hos8CRNF7SdyU9KWmBpG9LGl+M4MxsaHCeMbN8y6UH5zZgNfAB4CxgDXB7IYPqLw8yNhu0Bk2eMbPBIZcCZ4+I+HJEvJK+vgKMLXRg/eFBxmaD1qDJM2Y2OOQyyPhBSWcDd6SfzwJ+WbiQzGwIcp7pg28YarZ7eixwJG0ieeidgM+SPAQPkl6fzYD3LjMbEOcZMyuUHguciBhdzEDMbOgZaJ6RdCNwOrA6Ig7tZv6JwC+AV9JJP4uIqweyTTMbHHK6D46k9wEnpB8fioh7CheSmQ1F/cwzc4FrgR/10uYPEXH6AMMzs0Eml8vEvwZcCjybvi5Np5mZ5UV/80xEPAysL3B4ZjYI5dKDcxowKyI6ACTdBPwJuKyQgfWHpDOAM6ZPn17qUKwLD5C0PhQyzxwv6SlgBfD5iFjcXaP00RAXAUydOjUPmzWzUsr1UQ1j2XGUVLbXYEfE3cDdDQ0NHyt1LJXCj5WwIipEnnkS2C8iNks6Dfg5cGB3DSNiDjAHoKGhIfK0/X7pbr8zs92TS4Hz78CfJD1IcqXDCcAXCxqVmQ01BckzEbEx6/29kr4vaUJErB3ous2svPVa4EgS8AgwGziGJPH8c0S8UYTYzGwIKGSekbQXsCoiQtKxJOMO1w10vWZW/notcNKk8POIOBq4q0gxmdkQMpA8I+lW4ERggqTlJPfNqUnXex3JDQM/IakN2AqcHRElPf002Pn0mQ0WuZyiekzSMRHxRMGjMbOhql95JiLO6WP+tSSXkZvZEJNLgfMu4GJJS4Emku7jiIjDCxmYDW0e3DzkOM+YWV7lUuCcWvAozJdR21DnPGNmedXbs6iGAxcD04FngBsioq1YgZlZ5XOeMbNC6a0H5yagFfgDydHVTJI7jVoJ+JRN3/wbDUrOM2ZWEL0VODMj4jAASTcAjxcnpP6rpDsZ+0oFGyIGXZ4xs8GhtwKntfNNRLQlt6oob76TsdmgM+jyTLnwuD2z3vVW4BwhqfMuoALq0s+dVzeMKXh0VnD5SJJOtDYAzjNmVhA9FjgRkSlmIGY29DjPmFmhVJU6ADMzM7N8y/Vp4jYIeaCymZkNVS5wzMyGCI+Xs6HEBc4g5mRVHP6dzcwGH4/BMbNBS9KNklZLWtTDfEn6jqQlkp6WdFSxYzSz0nCBY2aD2VzglF7mnwocmL4uAn5QhJjMrAy4wDGzQSsiHgbW99LkTOBHkXgMGCtp7+JEZ2al5DE4ZlbJ9gWWZX1enk5b2bWhpItIenmYOnVqzhvw1Ypm5ckFjplVsu6e/RDdNYyIOcAcgIaGhm7bDCYuvGyoK/sCR9IBwOVAfUScVep4zGxQWQ5Myfo8GVhRolgGJRdKNlgVtMCRdCNwOrA6Ig7Nmn4K8G0gA/wwIr7W0zoi4mXgAkl3FjJWM6tIdwGXSLoNOA5ojIhdTk9VAhciZjsrdA/OXOBa4EedEyRlgO8B7yY5unpC0l0kxc6/d1n+7yNidYFjrBidCc6JzoYKSbcCJwITJC0HrgBqACLiOuBe4DRgCbAF+LvSRGpmxVbQAiciHpY0rcvkY4Elac8M6ZHVmRHx7yS9Pf3S3wGCZjZ4RcQ5fcwP4B+KFI6ZlZFSXCbe01UN3ZI0XtJ1wJGSvthTu4iYExENEdEwceLE/EVrZmZmg04pBhnnfFUDQESsAy4uXDhmZmZWaUrRg1OwqxoknSFpTmNjYz5WZ2ZmZoNUKXpwngAOlLQ/8DpwNvDhfKw4Iu4G7m5oaPhYPtZnZjaU+EIFqyQF7cFJr3B4FJghabmkCyKiDbgEuB94DrgjIhYXMg4zMzMbWgp9FVW3VzhExL0kl2/mlaQzgDOmT5+e71WbmZnZIFJRD9uMiLsj4qL6+vpSh2JmZmYlVFEFjpmZmRlUWIHjq6jMzMwMBsHDNneHr6IyM8udr5aySlZRPThmZmZmUGEFjk9RmZmZGVRYgeOrqMyGHkmnSHpe0hJJl3Uz/3xJayQtTF8XliJOMyuuihqDY2ZDi6QM8D3g3SSPgXlC0l0R8WyXprdHxCVFD9DMSqaienDMbMg5FlgSES9HRAtwG3BmiWMyszLgAsfMBrN9gWVZn5en07r6gKSnJd0paUo385F0kaT5kuavWbOmELGaWRFVVIHjQcZmQ466mRZdPt8NTIuIw4HfADd1t6KImBMRDRHRMHHixDyHaWbFVlEFjgcZmw05y4HsHpnJwIrsBhGxLiKa04/XA0cXKTYzK6GKKnDMbMh5AjhQ0v6SaoGzgbuyG0jaO+vj+4DnihifmZWIr6Iys0ErItokXQLcD2SAGyNisaSrgfkRcRfwaUnvA9qA9cD5JQvYzIrGBY6ZDWoRcS9wb5dpX8p6/0Xgi8WOy8xKq6JOUXmQsZmZmUGFFTgeZGxmZmZQYQWOmZmZGbjAMTMzswrkAsfMzMwqjgscMzMzqzgucMzMzKziVNR9cCSdAZwxffr0UodiNmBXXXXVTp+vuOKKEkWyq3KOzcwMKqwHx5eJm5mZGVRYgWNmZmYGLnDMzMysArnAMTMzs4rjAsfMBjVJp0h6XtISSZd1M3+YpNvT+fMkTSt+lGZWbC5wzGzQkpQBvgecCswEzpE0s0uzC4A3I2I68E3g68WN0sxKwQWOmQ1mxwJLIuLliGgBbgPO7NLmTOCm9P2dwF9JUhFjNLMSUESUOoa8k7QGeDWHphOAtQUOJx8GS5wweGJ1nPm1O3HuFxET87FRSWcBp0TEhenn84DjIuKSrDaL0jbL088vpW3WdlnXRcBF6ccZwPM5hlGJf0el5Djza7DECbnHmlMOqagb/XXKNXlKmh8RDYWOZ6AGS5wweGJ1nPlVwji764npetSWSxsiYg4wZ7cD8N9RXjnO/BoscUL+Y/UpKjMbzJYDU7I+TwZW9NRGUjVQD6wvSnRmVjIucMxsMHsCOFDS/pJqgbOBu7q0uQv4aPr+LOB3UYnn5s1sJxV5imo37HZ3dIkMljhh8MTqOPOrJHFGRJukS4D7gQxwY0QslnQ1MD8i7gJuAH4saQlJz83ZeQ7Df0f55Tjza7DECXmOtSIHGZuZmdnQ5lNUZmZmVnFc4JiZmVnFGRIFzmC5lXsOcX5W0rOSnpb0W0n7lWOcWe3OkhSSSnKJYi5xSvqb9DddLOmWYseYFUdff/dTJT0o6U/p3/9pJYjxRkmr0/vKdDdfkr6TfoenJR1V7BgLyXmkuHFmtXMeycFgyCFpHMXLIxFR0S+SgYcvAQcAtcBTwMwubT4JXJe+Pxu4vUzjfBcwIn3/iXKNM203GngYeAxoKMc4gQOBPwHj0s+Tih3nbsQ6B/hE+n4msLQEcZ4AHAUs6mH+acB9JPedmQ3MK8XvWcK/I+eRPMaZtnMeyV+cJc8h6baLlkeGQg/OYLmVe59xRsSDEbEl/fgYyT0/ii2X3xPgy8B/ANuKGVyWXOL8GPC9iHgTICJWFznGTrnEGsCY9H09u97rpeAi4mF6v3/MmcCPIvEYMFbS3sWJruCcR/LLeSS/BkUOgeLmkaFQ4OwLLMv6vDyd1m2biGgDGoHxRYmumxhS3cWZ7QKSKrfY+oxT0pHAlIi4p5iBdZHL73kQcJCkP0p6TNIpRYtuZ7nEeiXwEUnLgXuBTxUntN2yu/+GBxPnkfxyHsmvSskhkMc8MhTug5O3W7kXWM4xSPoI0AC8s6ARda/XOCVVkTyx+fxiBdSDXH7PapLu5RNJjmL/IOnQiNhQ4Ni6yiXWc4C5EfENSceT3Nfl0IjoKHx4OSuH/ahQnEfyy3kkvyolh0Ae96Oh0IMzWG7lnkucSDoJuBx4X0Q0Fym2bH3FORo4FHhI0lKSc6h3lWCAYK5/77+IiNaIeIXk4YoHFim+rnH0FesFwB0AEfEoMJzkwXTlJKd/w4OU80h+OY/kV6XkEMhnHinFIKNivkjwltpZAAAgAElEQVSq65eB/dkx+OqQLm3+gZ0HB95RpnEeSTKQ7MBy/j27tH+I0gwOzOX3PAW4KX0/gaRbdHyZxnofcH76/uB0h1cJYp1Gz4MD38vOgwMfL3Z8Jf47ch7JY5xd2juPDDzOssgh6faLkkeK/sVK9GOeBryQ7tSXp9OuJjl6gaSS/R9gCfA4cECZxvkbYBWwMH3dVY5xdmlbksSU4+8p4BrgWeAZ4Owy/jc6E/hjmrgWAu8pQYy3AiuBVpKjrAuAi4GLs37P76Xf4ZlS/b2X8O/IeSSPcXZp6zwy8DhLnkPSOIqWR/yoBjMzM6s4Q2EMjpmZmQ0xLnDMzMys4rjAMTMzs4rjAsfMzMwqjgscMzMzqzgucKxHktolLUyfkPtU+hTifv+bkbS5h+lzJZ3Vx7LnS9qnv9s2s+JzDrFSGgqParD+2xoRswAkTQJuIbk76xUliOV8YBFFuDOupExEtBd6O2ZDgHOIlYx7cCwnkTwh9yLgEiWGS/pvSc9I+pOkd8H2o6RrO5eTdI+kE7M+f0PSk5J+K2li1+1IOlrS7yUtkHS/pL3TI7MG4Ob0aLCuyzKflvSspKcl3ZZOG5UV39OSPpBOPyedtkjS17PWsVnS1ZLmAcd3F0c+f0+zocY5xDmk2FzgWM4i4mWSfzOTSG5LT0QcRvIQt5skDe9jFSOBJyPiKOD3dDmKk1QDfBc4KyKOBm4EvhoRdwLzgXMjYlZEbO2y3suAIyPicJI7YgL8G9AYEYel03+Xdk9/HfhLYBZwjKS/zoptUUQcB8zrLo7cfiUz64lziBWTT1HZ7up80uvbSXZgIuLPkl4FDupj2Q7g9vT9T4CfdZk/g+QBew9IAsiQ3NK7L0+THJn9HPh5Ou0kkucBkcb4pqQTgIciYg2ApJuBE9Jl2oGfDjAOM+ubc4gVhQscy5mkA0h24tV0/0h7gDZ27hns7Yis63NCBCyOiON3M7T3kiSZ9wH/JumQdF3drb8n27LOmfc3DjPrhXOIFZNPUVlO0nPd1wHXRvIAs4eBc9N5BwFTgeeBpcAsSVWSpgDHZq2mCui80uHDwCNdNvM8MFHS8el6a9JEA7AJGN1NXFXAlIh4EPgnYCwwCvg1cElWu3Ek3cbvlDRBUoakW/z33Xzd3uIws35wDnEOKTb34Fhv6iQtBGpIjqp+TPLUXIDvA9dJeiadd35ENEv6I/AKyVNgFwFPZq2vCThE0gKgEfhQ9sYioiUdDPgdSfUk/z6/BSwG5qbb2wocn3UOPQP8JG0v4JsRsUHSV4DvSVpEcsR4VUT8TNIXgQfTtvdGxC+6fuk+4jCz3DmHOIeUjJ8mbmZmZhXHp6jMzMys4rjAMTMzs4rjAsfMzMwqjgscMzMzqzgucMzMzKziuMAxMzOziuMCx8zMzCqOCxwzMzOrOC5wzMzMrOK4wDEzM7OK4wLHzMzMKo4LHDMzM6s4LnCsoCQ9JOnCHuZNkxSS/FR7swrWdV+XdJ+kjxZhu1dK+slutA9J0/OwXee9MuACZ4iTtFTSVkmbJG2Q9H+SLpZUVv82JJ0v6ZFSx2FWqbJywWZJqyT9t6RRhdhWRJwaETflGNNJhYhhMHDeG5iy+k/MSuaMiBgN7Ad8Dfhn4IbShlR+JGVKHYNZgZ0REaOAo4BjgH/t2kAJ/99hZc//SG27iGiMiLuADwEflXQogKR6ST+StEbSq5L+tTPBde0C7qH79S2SHpfUKOkXkvbobvvpdm6QtFLS65K+Iikj6WDgOuD49OhyQw/Lny/p5bQ36hVJ52bN+5ik59J5z0o6Kp1+cNqdvEHSYknvy1pmrqQfSLpXUhPwLknDJP2npNfSo9zrJNX19zc3K0cR8TpwH9CZAx6S9FVJfwS2AAf0tL+m7TPpfrJW0svAe7PX3/UUTnf7p6QfA1OBu9P9/p/StrPTnuYNkp6SdGLWevaX9Pt0PQ8AE3r7npK+kMa/QtLfd5lXlnlP0mnpb7Qpbf/53r7jUOYCx3YREY8Dy4F3pJO+C9QDBwDvBP4W+LvdWOXfAn8P7AO0Ad/pod1N6fzpwJHAe4ALI+I54GLg0YgYFRFjuy4oaWS63lPT3qi/ABam8z4IXJnGMQZ4H7BOUg1wN/BrYBLwKeBmSTOyVv1h4KvAaOAR4OvAQcCsNM59gS/txm9hVvYkTQFOA/6UNfk84CKSfeFVethf07YfA05PpzcAZ/WyrW73z4g4D3iNtFcpIv5D0r7AL4GvAHsAnwd+KmliurpbgAUkhc2XgR7H+Ug6JV3+3cCBQNdTYeWa924APp7muUOB3+1GTENLRPg1hF/AUuCkbqY/BlwOZIBmYGbWvI8DD6XvrwR+kjVvGhBAdfr5IeBrWfNnAi3pere3BfZMt1OX1fYc4MH0/fnAI718j5HABuAD2etI590PXNrNMu8A3gCqsqbdClyZvp8L/ChrnoAm4C1Z044HXin136Nffg30leaCzel+9Crw/c59Kd2Pr85q29f++jvg4qx57+kmL1yYvu92/8yK6aSsz/8M/LhLm/tJCpmpJIXCyKx5t2Tnpy7L3dglNx2Uxji9nPMeSdH3cWBMqf/NlPvLo7itJ/sC60mOhGpJEl6nV9P5uVrWZdkadu063i+dvlJS57SqLsv2KCKaJH2I5IjshrQr/XMR8WdgCvBSN4vtAyyLiI4u8WV/t+ztTwRGAAuyYhRJ0jKrBH8dEb/pYV72vtDX/roPu+73Pelp/+zOfsAHJZ2RNa0GeDDd5psR0dRlu1N6WNc+JL093cVYznnvAyRjo74m6Wngsoh4dDfiGjJc4NguJB1DsiM/AqwFWkl2xGfTJlOB19P3TST/6Xfaq5tVZieYqen61naZvozkSGZCRLR1s47oK+6IuB+4Px0T8xXgepJemmXAW7pZZAUwRVJVVpEzFXihh+2uBbYCh0QyRsFsKMneF/raX1ey637fk572z67b7Gz744j4WNeGkvYDxkkamVXkTO1mHbnEWLZ5LyKeAM5MT7FfAtxBz0XckOYxOLadpDGSTgduI+l+fSYi2kl2oK9KGp0mkc8CnQPsFgInSJoqqR74Yjer/oikmZJGAFcDd6br3S4iVpKMhflGGkeVpLdIemfaZBUwWVJtD7HvKel96VicZpKu9s5t/BD4vKSjlZiefo95JInqnyTVpIMVz0i//y7SIuh64JuSJqXb3VfSyT3+qGYVKIf99Q7g05ImSxoHXNbL6nraPyHZ7w/IavsT4AxJJ6cDcYdLOlHS5Ih4FZgPXCWpVtLbSfbnntwBnJ+Vm67I+n5lmffS73WupPqIaAU2siPPWRcucAySqxQ2kRxNXA5cw86D6T5FUgi8TNKrcwvJ+Wsi4gHgduBpku7ee7pZ/49JxrO8AQwHPt1DHH9L0i38LPAmcCewdzrvd8Bi4A1Ja7tZtgr4HEmvzHqSQYGfTGP8H5KBwrcAm4CfA3tERAvJgMZTSY6svg/8bXpaqyf/DCwBHpO0EfgNMKOX9maVqrf99XqSsTFPAU8CP+tpJT3tn+nsfwf+VckVU5+PiGXAmcC/AGtIctYX2PF/2YeB40hywBXAj3rZ7n3At0hyyxJ2HaxbrnnvPGBpmn8uBj7S03cc6pQOWjIzMzOrGO7BMTMzs4rjAsfMzMwqjgscMzMzqzgucMzMzKziVOR9cCZMmBDTpk0rdRhm1o0FCxasjYiJfbcsHecQs/KVaw6pyAJn2rRpzJ8/v9RhmFk3JPV2V9uy4BxiVr5yzSE+RWVmZmYVxwWOmZmZVRwXOGZmZlZxKnIMjplZvrW2trJ8+XK2bdtW6lDK2vDhw5k8eTI1NTWlDsWGOBc4VrZe/+2d0NHOuL94LyNGjip1ODbELV++nNGjRzNt2jQklTqcshQRrFu3juXLl7P//vuXOhweeeQRNmzYAMDpp59e4mis2HyKyspXWwu0t/HmU/9X6kjM2LZtG+PHj3dx0wtJjB8/vmx6uTqLGxuaXOBY+RoxGqpr2LvhnaWOxAzAxU0Oyuk3OuGEEwCora0tcSRWCj5FZWVr33ecUeoQzGwQGzNmjE9NDWHuwTEzK4CIYNWqVbzwwgusWrWKiMjLer/61a9yyCGHcPjhhzNr1izmzZvHhRdeyLPPPtuv9S1dupRDDz20zza33HJLv9ZvViruwTEzy7OI4LHHHmPDhg20t7eTyWQYO3Yss2fPHtApnEcffZR77rmHJ598kmHDhrF27VpaWlr44Q9/mMfod9VZ4Hz4wx8u6HbM8sk9OGZmebZ69ertxQ1Ae3s7GzZsYPXq1QNa78qVK5kwYQLDhg0DYMKECeyzzz6ceOKJ2x8tMWrUKC6//HKOOOIIZs+ezapVqwB46aWXmD17Nscccwxf+tKXGDVq1ysT29vb+cIXvsAxxxzD4Ycfzn/9138BcNlll/GHP/yBWbNm8c1vfnNA38GsWFzgmJnlWWNj4/biplN7ezsbN24c0Hrf8573sGzZMg466CA++clP8vvf/36XNk1NTcyePZunnnqKE044geuvvx6ASy+9lEsvvZQnnniCffbZp9v133DDDdTX1/PEE0/wxBNPcP311/PKK6/wta99jXe84x0sXLiQf/zHfxzQdzArFhc4ZmZ5Vl9fTyaT2WlaJpNhzJgxA1rvqFGjWLBgAXPmzGHixIl86EMfYu7cuTu1qa2t3T6w9uijj2bp0qVAcnrrgx/8IECPp5p+/etf86Mf/YhZs2Zx3HHHsW7dOl588cUBxWxWKh6DY2aWZ5MmTWLs2LG7jMGZNGnSgNedyWQ48cQTOfHEEznssMO46aabdppfU1OzfZxPJpOhra0t53VHBN/97nc5+eSTd5r+0EMPDThus2JzD46ZWZ5JYvbs2Rx11FHMmDGDo446asADjAGef/75nXpUFi5cyH777ZfTsrNnz+anP/0pALfddlu3bU4++WR+8IMf0NraCsALL7xAU1MTo0ePZtOmTQOK3azYXOCYmRWAJPbcc08OPPBA9txzz7zcAG/z5s189KMfZebMmRx++OE8++yzXHnllTkt+61vfYtrrrmGY489lpUrV1JfX79LmwsvvJCZM2dy1FFHceihh/Lxj3+ctrY2Dj/8cKqrqzniiCM8yNgGDeXr3gzlpKGhITqvKDCz8iJpQUQ0lDqO3nSXQ5577jkOPvjgEkU0cFu2bKGurg5J3Hbbbdx666384he/KMi2BvtvZeUt1xziMThmZkPAggULuOSSS4gIxo4dy4033ljqkMwKygWOmdkQ8I53vIOnnnqq1GGYFY3H4JiZmVnFcYFjZmZmFccFjpmZmVUcFzhmZmZWcVzgmJkVQHt7O/fccw9f/vKXueeee3Z5NlV/ZDIZZs2axSGHHMIRRxzBNddcQ0dHR7/X190DNwHOP/987rzzzl6XnTt3LitWrOj3ts0KzVdRmZWhzRvWE4i6kSPJZKppa2mmZnjd9vntba0gQQSZ6hoigs0b1tHe0kLN8Doy1dVsWr+GmmF1tLU0M2LMWEaMGUvLtq1sa9rI6D0m5eXGc9a99vZ2Tj75ZObNm0dTUxMjR47kuOOO4/7779/lGVW7o66ujoULFwLJE8s//OEP09jYyFVXXZWv0HM2d+5cDj300B4f3GmltXXrVl5++WWmTp3K8OHDaWlpobq6evuT6COCLVu2UFtbiySqq6vZsmULL774ImPGjGHMmDGsXr2azZs3M3LkSLZs2cIRRxxBdXU1r776KiNHjmTixIkl/pa9c4FjVmIRQURABK2tLax77WXoaIOODjapCgR0tIPSDtdMBtrbkwJHYsTYCVRlqti8eiVEezKvvR2qqmitSpZp2bCWDW2tEAE1tWxe9TqZEWOYOPUAqqrckZtv9913H/PmzWPz5s1AcgfiefPmcd99921/EOZATZo0iTlz5nDMMcdw5ZVX0tzczCc+8Qnmz59PdXU111xzDe9617uYO3cu8+fP59prrwXg9NNP5/Of/zwnnngiAJ/73Od48MEHGTduHLfddtsu/2ktWLCAz372s2zevJkJEyYwd+5c/vjHPzJ//nzOPfdc6urqePTRR6mrq+saohVRW1sbmUyG1tZWXn75ZZYsWQLAK6+8slM7SWQyGUaOHEljYyMAw4cP5y1veQuvvfZar4/keOONN3aZ9va3v52xY8fm8ZvkjwscsxLpaG+ntWUba5e+CI1roXY4qqlNCpkIaE+eB0RVZntvDe0d0NGaFETVw1B1hi1rVkB7G7S1psVPGxDQAaBk+UwmKX6SDYNE+5ZNvPH80+x54CFkqmtK9CtUpj/96U80NTXtNK2pqYmFCxfmrcABOOCAA+jo6GD16tX85Cc/AeCZZ57hz3/+M+95z3t44YUXel2+qamJo446im984xtcffXVXHXVVdsLIYDW1lY+9alP8Ytf/IKJEydy++23c/nll3PjjTdy7bXX8p//+Z80NJT1TakrWkTQ2trKggULWLduXc7LtLW1bS9uALZt28bixYv7FcMjjzzC+PHjOf744/u1fCG5wDErsjeWLqGjaWPyIZOBbU2wbXNS0NSMh+iAzieopD07EEkB09FOtLcRq16D2joYNxHV1hHNW5Nen9rhSFXQ0ZEsJ4EyQPpndOxYrwIQq5Y8y7jJB1A7vI6qTManrvLgyCOPZOTIkdt7cABGjhzJrFmz8r6tzsftPPLII3zqU58C4K1vfSv77bdfnwVOVVUVH/rQhwD4yEc+wvvf//6d5j///PMsWrSId7/73UBy6m3vvffO91ew3dTW1savfvWrUoex3bp163jggQc46aSTaG5uZvjw4aUOCXCBY1Z0HVubdvTOtG5DATFiDGQ6e1EEBNHSnLQbNgJVZY3bUBXUDEt6Yt5cRYwev6NwaWsjVAXNW6B2GMpUp6tLi52qaiDtHYrM9m2++dqStBiqgupqJk2dTnVtbZF+kcpz6qmnctxxx+0yBufUU0/N63ZefvllMpkMkyZNoqfnClZXV+80EHnbtm09rq9rcRsRHHLIITz66KP5CdjyIh8D1vOtubmZX/7yl0Dyb27ixIkcffTRJY3JJ9/NimzS9JmAkgKlrRU6WlF1DcpkIFMNVem4m/a29NRTC7Q1J8uoClXXUrXnfmjsRBg2Empqkp6gqgyqqkradrRtPxW1/ZRVdOzcO9TRnhQ+He3J9I50flsra5cvLdXPUxEymQz3338/t956K1dffTW33nrrgAcYd7VmzRouvvhiLrnkEiRxwgkncPPNNwPwwgsv8NprrzFjxgymTZvGwoUL6ejoYNmyZTz++OPb19HR0bH9aqlbbrmFt7/97TttY8aMGaxZs2Z7gdPa2rr9VMbo0aN7Ha9hhTNs2LBunwZfLtra2li5ciVtbW0ljcM9OGZFtm1zIwyrg84emiDpOeksRjo6kgKlbiQg1FncREfyZ3qUrbpRRFUGtm1JemJatxGZ9JilZhhU16brZkchIyVFlKqSV1V6Ooudj/7HlPnVEYNBJpPh9NNPz+uYm61btzJr1ixaW1uprq7mvPPO47Of/SwAn/zkJ7n44os57LDDqK6uZu7cuQwbNoy3ve1t7L///hx22GEceuihHHXUUdvXN3LkSBYvXszRRx9NfX09t99++07bq62t5c477+TTn/40jY2NtLW18ZnPfIZDDjmE888/n4svvtiDjEukra0NST323JVadXU11dWlLTFUrj/OQDQ0NMT8+fNLHYZZt9Yte5nmpk2M3GMio/eYxBvPLkiKj+zBxB1toOrkc0faHd2cDFqNqurkFFR0QHsb0bJte/GjkfXJetLKRtU1yfic9jYYPRYpHXBcXUNS+bDj6ixg4vSZ1NQOK+j3l7QgIsp6ZGp3OeS5557j4IMPLlFEg4t/q8Lq6Ojg3nvvRRLvfe97WbRoEUuXLi11WNvls6jvTq45xD04ZkU2du8pbN24gbox49jwenoJZ1UmLTQi6VWpGpa87+xdkYjqWmjZBi1N0JqOo6iqTnpktjUlg46ra5Jl2lqTPpmqqqSXqKM96TGqG0VyeqxjR69RJ4nG1SuZMHlaEX8NM9tdVVVVNDQ0kMlkWLt2bcFuuNjR0UFE7Pap1fb29ryeju0vj8ExK7JMdQ2j9phIS9NmtjWuzzptpB2njiQ6x9xQVU20t0FTI7Q2J70xNbVJsTKyHo3dE9VPhFHj0uKmhWjaCC8sJNavSgcYVyVFVNrrk1ydlQ48ldJxP2LEmPK8n4WZ7WyvvfZi3LhxvPLKK7S0tBRkG//2b//GBRdcsFvrr66uLoviBtyDY1Z0m1a9zqbVr2eNfanacaqpKoOqa3cMEs7UbL+h3/YCqK0FOgQ1w6G9NamFho+A5q1JAVRdu2MMTnVNUgxVZdCwuqxtpgVUkNwUMJOhetQYFzh9iAhfRt+HShz2UG5aWlr4v//7v51uQ1AIe+yxB1u3bt2tm4G+853vLGBEu8cFjlmRtGzbwqbVK2jb2rSj9yQVHR3QsjW9yV/7jhqkowMymeTeNqP3AERsqUpmZtLTUZB8rq5J1jtsJJkxe1B/5DvZsPyl5KqqqkzWaTCScTiw4z47UUX9+D2L80MMUsOHD2fdunWMHz/eRU4PIoJ169aVzX1QKtGLL77I+vXrC17cZDIZ/t//+3/U1dXtcjfknowePbqsBpu7wDErgpZtW1n7wqLk9JCgdo9JDBs5BiIYMXYPVFXFykULkhv+dZ5S6jpGBoBAo8alY2q2kozPSYoVZarT+9ok99Cprqtj4owjaN22JRnrs/2mgZ2DiztPg4kJBxxMrf9T6tXkyZNZvnw5a9asKXUoZW348OFMnjy51GFUpMcff5zVq1cDyRVu9fX1TJ06lY6ODvbdd1+WLl3KokWL8rKt9vZ22tramDlzJuPGjWPFihXdPqqhU21tbVn13oALnAHbunUL6x67H4aNYPJfnFzqcKxMrX/l+fRqqOQ+NJnaYQwbOZra9AGarVs2w5bGZKBw17ExkPbsdGwvZoBkgHG0s71gSQcjA1SPGkumppZMTS01w+toXPU60day81ic6ICqaqqGj6C6xo9q6EtNTQ37779/Qdb9q1/9ira2Nvbee++S3xzNytPq1au3FzeQnKaaOHEikyZN2j7mZe3atXnbXm1tLXvssQeS2Geffairq+u1wJk6dWretp0vLnCydLS10d7eTs2w3C+TXbfoCdiyMbm6xSwVHR20tTSzec0KtjWuT05BxY7TSVvXr6Vl61b23P8gADa//nLy7wjBiDHQunV7ARLJChE7LiOP9LlSqqpOrpKqHp6Mu4mAqgxtW5tobtpIXf14mtavoXbESJo3d97oL2sMDjBh8jSqymRQYCVobGxk5MiROd8D5M0339x+Q7SVK1cWMjQbZFpaWmhsbOSpp56itps7iz/77LMMHz58+xPdV61alZftSqKlpYVXX32Vgw46iI0bN/LGG28wbNgwmpubd2k/evRo3vrWt+Zl2/nkAifV2tzMmuXJecYxE/ZkVP04OjraqarqPvG3t7WybsUyqupG0lE/EYaPKma4VuZWvfI8Hdu2JHcqbu98hlQ7sWQR2m8GGjGaUXtM2N5+1NQD2bJhHdTWEW3NoAyqGZ4UMmuWwYh6GD2O5JRSJL09WzYSo8Yldy/efiZrx0M6tzZuoAPYuHLZjvvsVGWd9lIyLqcq4zSQL7/97W/ZunUrkNwLpL29HUk9DtJ89dVXd3rIoXtvrNPKlStZsGDB9s+dj9hYuHAhL730Eu9///sZPXo0e+65Y+zcXnvtlZciuXOgeHNzM9u2bePJJ5/sdcxPqW/o15PyjKrIWluaaVybXk4LbFyzko1rkn8kNcNGMnbSntQM23l8wvo3ltO2uRGig6pxe7HnjMOKHreVl4gOtry5jmGjxpCpqSVam6mqqaW9cR3RvIVYsRReeJrYtgVqa9jS1kLTsiVUZzK0tzQzbNIUWrY2EavXJMVI/fCkCKlOH8UQQFXa81I7AlpbkqunqjLp6aq0Z0YCQXPTJpq3NkF0PqoBtt8ZIh3fM2Hagbt1hYR1LyJYvnz59uIG4N57793+DKgZM2aw//777/QfwcaNG3nmmWe2f95zzz39IEujsbGRTZs2MWLECGpqaqiqqto+HgaSR2o899xzvO1tb6Ouro7f/va3TJkyhZUrV7LvvvtSV1e307/DgfrNb37T6/wDDzyQGTNm5G17+TTkC5y1r79Gy7YtO57Y3Cl927qtiTXLXk7/c8mQqa5GgrYtm5M7y0pQU8uqpS8CMGxUPeP32rfo38NKq3XbVta/+iLtba0MHzmaqtZtRNNG2iOSp39vWA011bD/QbDnVKJ5Ky2qgpZttHYk42j2Oew4VJVhTXsLLVubINqTS8bHT056ZdpbydTU/f/t3Xl0W/d54P3vvRe42AmQBHeRlExK1mbJ1mJbjhNvtZLYipM2ruOmzaTTpD1vp30zM017ZtJ523PSt2ne9k1z2nnTc9q0qTvpyTR2Js3JWEnqJfJuS7YWa99IiSIpLgIJkiCx4977/nEBkDQpiZK4As/nHB4RF8DFDzQMPPj9nt/zEGppJ5tKEBvss7d4FxvvFRpqapOzNZDP28nP4IAd3GhOGtrXy26geRCLxXj99ddnHJ/a4PLs2bOcPXsWAL/fj6qqM/r0DA4OsnfvXgDuu+8+qqqqFnDUYjk6e/YsFy9eJJfL0d7eDjBjSeh3fud36O7uRtf14qzOhQsXsCyLiYkJHnnkEeLxOK+88sqsj+F2uwkGg9xxxx0cPnyYaDR60+PdunUrzc3NN33/hVb2Ac6M4KYY41hTtuCS7xNkYhSqwlrWZE8f1VG8X3pijGS8ArfXLx8eZSI+OszY5S67Bg2QGp4ozgZaRs5uq5BJ2/kv1XUw1AujOlbDGhTV3t4daG4vdgx3VFSRzeWL8aWTxZ5RimVhxGOMXOqgorF1MmlZU+wgptBUU7ugLcQAACAASURBVHNMJhIDaLp9O1UjvHotukt2S82njo6OG7r9XLb3Hj16lPvvvx+HwyHvI2XAsiz279/P8PBw8djVXletra3s3buXb3zjG/zVX/0Va9aswbIsfD4f69bZOX0+n++qfapSqRTpdJqenp5pQfhcNTU1cdddd93w/ZZC2Qc4M2Zu7INTttTmvxVbTOYuqA77uOafrAhbeCEpCiOXLxXrjTi9PmqaWhfnuYhFYVkWka7z5BITdq+nTHIymJi6DduysLIZO8BxeSY7g6v5JpeJMahtBRV7pgZIxydIpzP4G1fjCQSJ9V3CNE00txdvsIroxTOYuSyj3R35hOJ8403dhT/URM7IkYyNoukualtuK77JDXV34tBdEtwsgFsJQKqqqmb9Bh2Px3nhhRcAuyz/vffeKzM6JSYWi/Hmm2+i6zqZTOaGgo21a9dy6dIlQqHJwpwej6e4BHrw4EH8fj/r1q3DNE06OztRVZW2tjYuXrxINBotzihO1draiq7rDAwMEI/Hueuuu4rLpv39/Zw4cYLW1pXzeVbWzTYNI8fgxXP2BTMf0Ch84MOKKbVI7H9d/grSyUQxmXMaRZ1xP90XICxBzopkmQbx/m68dc2oDgdjV/qJX7lsz9YUqgxPfb2Yhn3ZNLDi45jvvAxVYRS3G2qb7OtT+fVxhxOqG2F00K40HKpB012Ybj+BhlZ8wUpUTbN3Yw0N4quqYehShx00WWZxG3n1mvW4PN5Zx58cHyMVnyBY27Bscm1WarPNqzl+/DiXLl26ofN7vV5M0ywuMczFI488sqyKqIm5GxkZIZ1OU19fTyaTYf/+/cRisTnd9yc/+Qnvv/8+1dXVtLS08NhjjxWvO3LkCD09PTzxxBPFY5qm4XQ6URSFXbt24fXa7w1dXV3kcjlcLhdHjx4FKH4BcrvdPPLII7MG64ZhcPz4cWpqamhqWh7pF9Jscw60wu4RiykfUle58ZSZnvTE2ORxRcvXIskzp/xu2R9+meEBIsk4Ne0b52voYpH0Hfg5TIyQGrlCRnPbr5NMMt9GwZFfnsx3/Fby1YWzaTvvpv8iRPrt2RqvBypCKL4K8PjtXlHD/RAdhMoa0CzIZjBMEyWXZaIrQcpXQc3tW0mORUmMRNB0Fw6Xm6xl4g4EqWpafc2xp+LjjPT3gKLgrwqjLnCX8HJVV1d3wwFOIpEo/q6q6py+vf/85z/nQx/6EJWVlTc8RrF0Tp06xYULFwB7l9zUnVEwOQN4tcmGgwcPcuTIEQCam5unBTh///d/T09PD6lUiqeeegqwA5LC7r19+/axbds2GhsbOXv2LLlcjl27dgF2Xaft27cTDodnPugU+/fvLwZoyyXAmauyDXBM0+BK94XZryzMwqjqZL7NVFZ+2co08+VEpkS9heUqhcmy+ChkUwmS42N4AsEFeDZivpmmwcCpI3YQo3tw164iM5SvMaF786+NfM6LmW9gmc+VsbdiayiVYXjkkxBeBeMjKHq+mJ7HZ/eRGonY28ir6uxzev24vX6MRAzF6Ub3BQAw8vk7DoeDysYWctkM7vx1V2NZJtHOk3YDzbpmnBLcLIje3l7ef//9WzrHjSxN7N+/n49+9KPLZjZOXNvRo0fp6ekB7ED28OHDM25zvVWUP/iDPyAWi5FMJvH5fNOue/rpp/nmN7+JYUx+sW5ubmZ8fBxVVclms3i9XnK5HNlsFrfbTVVVFTt37qSiouK6M4KnT59mZGQEgC1btszpOS8nZRvgjEeHMXP5JSar0IAQ7AaE+dmawguvkI9TaHgI9oeaZdmzNzPebKwp/QwVO//CNIhePIvuD1Bz24aFfnriFg1eOAdYdoNKt4/YSHRydk51UrW6nWh+5xxGurhkZEF+Z50Lws2oar5TuF6b34mn2ktTioJ1+zb7MfIvrZq2TTjcHoYunUfVHARX3QaAt6oGRVFw+YOoDgeOOQQro70X7W3kikpFde0C/IWEYRjFqf7FfMyf/vSnsly1AhiGUQxuwE78HR8fL17etGkT586dI5udJdVhCq/XW1xmmsrhcPCRj3yEXbt24cxXIvf5fGzdupXLly9z6tQptm7dWszTWb9+fbFH2NTaOdfS2dkJ2DuvVuLrrWy/BsTH8ol9U5enCgcKszBW/tt58bJpXzbN6cHO1Ai80N8HJoOgwuWhXjLnj5CemNvaq1g63kA+ea+w3drI5l8PFuRSRM8ft2dfsEB3g6ZjxoaxLndAchzMLIqVn9nJBxpo+UrEWCiajqI5UBL2G56ruh6n18dQdyfZZJxcJkX/qcP0nTiIoqgEG1tRb6CYlqeyBjQHzkBICvktkPHx8SXrnH21LcBi+fjgLNvU4Abg5MmTswY33d3dxVmTaymUGYhEIrz00kuAvW07Go1y5MgR0uk0XV1d/OQnP2H//v20t7ffcI+wYNBecbjnnntu6H7LRfm+8ynq5DbwQp8fmD5rY5rAlG7NU2KVydsZdg5G/kPEFQrjD1SQGB0hOTY8eRtVBY8fNI2x0RFq/RUL/ATFrQjWNeCpCDLUeXKy/1Oh03fhdWPlACeoCpbmgKE+iMew3H4UT8CuT1PIzzKyWJmUfcxXYW8JT2WwBrvBHyLtD9J38iCKy4fuD+INVjF6yU6At6y5L2EUuP0VNN553zz8JcTVaAvU3kLTNO677z6uXLlCT0/PtHydAtM0sSxLtpAvY4qisGfPHvbv3z/nHlFjY2N86Utfor29nW984xtzus93v/td3nrrrWlFIisqKlizZg0dHR1YljXra2guPvzhD9/U/ZaLsg1wXB4P6Yl8RK0odpAyrdDfbNvH88eLwdHM3Jz0cD+ZsWHCq28nnZzAzGbsq+JjWMkJlLpWconxmecVy47u8dpJw8VyASqohdwsCxxOVN2Foih2VetgtT3j4w/a1+tuO6hJxOzZm0zKXubSdbu1hzeAUtti3z7PyqYJt9tLmOqa29GcLpzulTc1XA4CgWvnQd0swzB44403aG5upr29nWPHjs16uxMnTnDHHVJBfbmrqam5ZoCj6zo+n4+RkRH8fj+/8Au/cENNXZ966ina2tqm9YKqrKykubmZuro6IpFIsVdVuSnbJaqq+lW4pnywTPbnmfKvouZ/NCanbqYsTRUvq5M/5LcWj40Qbl2L7g8CCtb4KIyPQHpxmnIOX75Ef+dpDCN3/RuLq/JVFdaq8//d1XwRvXz7BDMVx0glQHOgVDeg3LYVJRWHiShm/0WsaH++KKQTgjUQqLJr3uR3ORCqhYlRrLFhUDSc/sm6Fu5ASIKbZe6BBx6Y1/MVlgQALl++TE1NDdu2bZu118/NFGm7EZlMhpdffnlaryxx49ra2mZtlFmgaVpxSUrTNH73d3+Xxx9/HLCXq/7sz/6Mrq6uq97/tttu45d/+ZenvUY2brR37Oq6TlNTU9nO9JXtDI6iKHgrKkjHx+zZm2I+TWHZ6oOBzBTFdffCrE++waHDicMdJJdOkRgdJjk+hpXL2J+NtasgVIu7pgHdPXvNkps1PjJEfGwUM19Jt1gBF4j2XaKmuW1eH6+cBOsacPr8jPZ0ovuCWJZBNhnHW12L2+MjevGM3VKhsEUcC0t32xWMRy/ZycSr1qH4Qii59OSMkJWvVJzL2Lk8WobK1nY819kdJZaXQCBARUXFnGuaXM/Y2BjV1dUMDw9jmiavvfbajJYOLS0tJJPJYin/+ZBMJjl9+jRXrlyZ8XgXL15k48aNZfshOR92797Nq6++Sjwe56GHHmLfvn2oqsquXbs4fvz4VXtHHT9+nHfeeYfNmzezevXq6z6O3+/ngQcekP9WeWUb4ADobg9uf5DUeL6ujaIA6vS6NlMVqxszZaZHndxdYxjkknF7CcMy8rGSlb+ZRt3GbYxHIzjyGe83qrDmbpkm8bERYqPDkBi3x+R0Td/NpWqARaC6/qYeS0zy+gN4N9w5+3WVNeSyGXK5HGYqDoDiDaKYBmbLersTeHwExe2b+bqyDHt1tLKOujvuxeGSrdwrUVtbG+fOnSMej8/L+aaW6/9gsLF69Wrq6+vp6ekp7py5UYX3kdHRUXp7e685OwD2rIJ8YN66Bx98sPj7nj17ir/fcccdnDlzBsMwGB8fxzCMYm2kj33sY6xevXra8tPV6Lo+7TFEmQc4muagqr6JvkLhPkuZrF8zW2LnbDsmFMXOrQBw6pO1UABPIEgyv1vL6fGTzaRIxEZIRq8wioWzopqqhlWz1rTIJOOgqIwNXCaXywImViaDZebs2RmXD0XTIJuyx+VwMrVTtL+xhYrKaxdwErcutGpyrdzIpEFVScXGGBvsRdU0LEW1g8+p+VqF14hllxMIb7hzWnBjZNLEI324/CGM5ASxjmPkdC81m3fi8voX8dmJuWhqaiIQCMzacPNWKYpCRUUFY2P2e1RLSwvHjh1jdHSUgYEBdF1nw4YNs+ZYZDIZkskkg4ODxeaMNzLTFAqFuPvuu6+5vCJuXWVlZbH4nmVZjI+P43a7OXjwINFolE2bNl33HLqus3v37mnHent7MQyDQCBAf38/PT09uFwuHnrooQV5HstRWQc4BZUNzfmKr0zWKjGZGeQoymQBwEJdnMJlyAdAZvGDKzkewxMKk4yNkE2niPZ12zczDMAiE48xcOEMTlXF4faABf5wHZlkgrGezsk8IMuEbBrj/DGIRlDWbUIJmOAL2gXj8onPTl8F4YZV8m1riWj5+jS+qjC6z0/k4lmUQL7qbC6b70OlTUnnUnEF/CQvX0Rv24yiqgy99W+k9r8Emgq3rQdUiA2DZRG53EnjE/8edYF274ibN99LVQWWZREKhUgkEuRyOd5+++3irI5hGCSTSQ4fPkxvby/ZbBaPx0N1dTWNjY0z2gGkUilSqRTf+9732LFjB1u3bp31MR9++OFZ666IhVcIaMHuKH/w4EEGBgaueR9VVamtreXKlSvU1tYSiUT4rd/6LZ5//nmefPJJfu3Xfo0DBw7wp3/6p6iqyosvvsgjjzyyGE9nyUmAA3h8ATztG+nrPJ3vyKzkg5gPFP0r5OUUfgrH9A8kghbq5eRMkrFCPYMp3+BdbgrRlGJZZFMJsqkE5HIkh/ux0vnmjQ4XZNPkThxh4rv/jKI7CWy5DeqbIbwKh8dH5ao1Nz1VLRaO0+Wm4fYtTEQjZFIp0iNX8hWOp752IDXYjZXLkOjthJErmOdPQWICgkHIZOwu5IOX7erHFdVYpmnvyBLLiqIofOQjHyGRSLBv3755PffUNhAfXLIquHLlCmD3POrr6+P48eMzbvNLv/RLxbor7733Hn/7t3+L0+mkra1tXvN5xPzZsWMHExMTdHV1MTg4OGuujmma9Pb2MjAwwE9/+lN+9KMf0dvbC8BLL73EZz/7Wf75n/8ZsNszlNPnhQQ4U3j8QZKjQ2i6jtPjs3NzCi0ZCiX5YeYHzAdnTOxytvbPlF1MFY0tZFNJ3G4fyfgYqfEYlqKAotlLT5FeCIXtb/oW4PQAFrnnf0hmLAlWkuo77qb2ic+jSKn2ZU9RFALVtWQzaSITY3atpELAnLPLB6B7IJfFHI9CPAY1teBug1AVuLzQdTofLGdpePxzaI7yeXNaidxudzF/YsOGDXR2dpLJZObt/F6vl7Vr15JIJAiFQhw7doxsNotpmiiKQm9vL0NDQ9NmZzRN4wc/+EExuHE4HHz729+e1tNILF9+v5/NmzfT399/zdvlcjk6Ozvp7+/n6aefRtd11q9fTyqVKlZU/trXvrbia9vciLLuJn49kd5LZFNx+0PJMCYbKqqqveSgajODm+I6F5Pf1At/43zirzsQIjUew65o68Dp9mD0nCN76FWUlvVU3LcbbyCE5nJjmSaXv/IFUidP4dv9OA1f+uNbfl5i8Y0M9JJKJbGS+UTUbAawsDJprDOHoNBOIZ2AhjW4a5txqCqJ/kuYg70QqCB420YCq6+fbLjclVo38WvJZrO88MIL8zCimdauXcv583a7kEAgQHV1NU8++SRnz55l79693HfffYRCIRRF4cyZM2zYsAFVVUmn07NuOxfLWzwe59ChQ0xMTEwrEaCqKq+99hr9/f08+eSTpFIpvF4vqqqyadMmuru7OXnyJIlEgh07drB9+/YZPa1WGukmPg+qG5sZH76C7vHh8dvbdzPJJENdZ/PtG/K7rmCyRcO0Vg35oKi45GWhOhx4K0IEqmsxcll0jw9VVTGqaxnPpvFu3IEentz5pKgqq/78mcV82mIBVNbbJdKz6RSRjlNYEyNQUQ2DPVjnT8D+IVjTjlJTj6b3U739QYZPvofZ1wXZBKRijKUmyGo6/qqaYiNOsbw5nU527txJJBLh9ttvLy4PzCW34mpCoRCWZdHe3k4oFMLj8RTzNr7+9a9z7NgxHnvssWm5eOvXr1+ythJifvh8Pj7ykY8AdpHHf/u3f6OiooLGxka+//3v093dzfPPP09dXR3RaJRvf/vbVFRU8Hd/93fFXVh9fX2kUilWr15NW1vbglXjXi5kBucmTIwMkY7H8IXriV44C4qCoruwpq6PK4o9y2NkAQXNGyDc0oZWRuufYibLNOn/4d/aW8qb2sFTgfXqj+H8aVjVjLJmLam39mMORPDceyeKy2nn37jc0LQGta4VVJWGrfet2GXKcprBuZpMJsPRo0epra1ldHSUnp4eFEXB5/MxMTEx631UVaW1tXVOu2pEaTt8+DDbt2/H4/Hw3HPP0d3dzZe//GVSqRQej4dkMomiKASDQUZHR4v327ZtG1/96lcBCIfD3HvvvUv1FG6JzOAsIH9lGH9lmJHLXaAoBBta8IWqsEyToctdZJNJVFXFGfBi5XIEG5qlIq0AYOxnz2L87Eew+U5cm+4hp3sxN92F5fOAvwKqG5g4cg4zk0NfP4rmcoLba3c194Ww4mPgCZBOTOCWfmYrlq7r7Ny5k1QqxcmTJ/H7/dx3333ouk5/fz8nT54klUrR1NTE6OgoDQ0NrFu3btaSEqK8WJbFxz72MVRVxe/3o2kaLS0tfOELX+Bv/uZvePrpp/nxj39MNBpldHQUl8tFOp3G4/EUZ4AAhoaGSr6fmfzfchMMI0difIx0ws7P0fKNNhVVpab5NnS3B0VVqVq1hvCadTcd3PT97F+48B8+zZH/81cWvCy7WHjjvZ3E4zHUhhaqH/8cjvrVdp2cdAoudEB0FOvYu3hXh/GuqkQL+O2dVKFqCIVR8tvQUVWiXeeW9smIWxaNRunt7cU0TXK5XLHeTENDA9u2bUNRFCorK3nooYdYv379TQU3haJxiqLwm7/5m/P9FMQiMwyDN998k82bN/PUU0/R1dWFYRicO3eOZ599FoDR0dFieYCNGzeSzWZZvXo1iqKwefPmaecr5HCVKpnBuQljV/pJTYwDCoHaBlz+AEYuR3J8DG9FiHDLbfMSGade+zkH9h7gZz1jPHrhEzz9P/8FZ1C+ta9ERnKCsQMvggbez3wR99pNxE68izI+hvnaC6QvX0EbjuFwa3g2toPPh7r5HkxfEGWoO9/7ymlXRHY4QbHzeZwu91I/NXETEokEb7/9NmA3Yyxs0x4cHERRFGpra2fk0dysQhrCP/zDP/Doo4/y1FNP3fI5xdI4e/YsY2Nj/N7v/R4f//jHuXz5Moqi8J3vfKfY0PNHP/pR8fa5XI5nnnmGn/zkJ+zbt29Gzs25c+dYt27doj6HxSQzODfBWxHCobvwhaoJhOtRFIWxoQFikX6S+bYP8/HGpKzfQpXuoNbtYPyNd9n/C/eT7uq85fOKxTdx9n0cbh+OQAh/cxsKdvkB89DrGIkksY4I4xcj0NRK9Rf+K56t92GpoDo0lIY2lMpae5v5lJ17ka7S/vZVytxuN42NjQSDQdavX091dTW5XI6DBw9SyP2Zj/eQWCyG2z0ZBH/mM5/h6aefxjCu0o5GLFuXLl3iL//yL5mYmKCxsRFN09A0jdOnT3P27Nni7dxuN1//+tf51re+xRe+8AWqqqr43Oc+xzPPPEM4PLO6fTQaXcynsahkBucmuH0BdI+Pwc5TxKNXcPsDpEaHQHPgCQSvf4I5qti2g3WP303z6Qv0dUVxjMc59yufpP25/42n+bZ5exyxsFKxUcYPvQqqRtPn/wu5dBLTMDj86V/D5dJov7sV37134rzjDpSNGxh97+cEtj+I1+Mjm04x0dsJugfFqU8pEjilB5pYcVRVZdu2bZw/f5433niDcDhMNBrFsqxZ2y7cLIfDwcMPP8xPf/rT4rFnn32WQ4cOlfzyRKn567/+a/7xH/+R2tpaPvnJT5LJZHjwwQfp6OgAwOPx8Oijj/Krv/qrvPXWW0QiEf74j/+YbDbL4cOHr1qPKZVKLebTWFQyg3OTTCNnV5W1TFJjI2CaaKo6b2X0UyMRkueOYfr9XL4UparWj9evozg1Bv7vL83LY4jFkRzqh5Z10Lqe2OWLRE6+R2JokHQ8QzqRBU3H/4kncN//MIwNYY4NM3ZwH1nDJLj6dnxrNkJyAk1357uP54q1lUpxF2Q5KSwrDA0NFfPsmpqa5u38X/3qVzlx4sSM4x0dHfzrv/7rvD2OWFiGYfChD32IL3/5y+zZs4fXXnuNd955Z1rxv/Xr1/PFL34Rr9fLM888w7e+9S3effdddF3ngQceYHx8fNqOqoLZjpUK2SZ+C/rOHIdM0q5xo2rUrds8b5VmxzuOcflfv0f/M/+LyJU4tS2VtKyv4dAbnQxnTX6j8xy6X2qhLFeWZTE+NEAi0o8RG0Zz+6hoXYuqqkwM9BJoWMXwj/6F1Olj6L/6m3hC1aRHhzGTMejpgGA1iqqgun24shkS7/4cpf0OlNbbJ3uUOd04glXUrlq91E/3hsg28UmxWGxak8477riD1tbWeTt/e3s7nZ2zL2urqipLVctcIpHg6NGjjI+Pk8lkaGhoYOPGjZw9exa3200wGOQTn/gEn/rUp9i+fTuVlZUMDQ3R0dGBYRjcfvvt5HI5enp6+Na3vsXFixf57ne/i98/vWnvRz/60RXVwmGu7yEyg3OTjGwGKzmOlRrHymZwV9XMaxn9QPsWzv5/3ydyJU5DlYe6KjcD71/mhb4Yb1+Z4MR//3/n7bHE/IrHRhg8f5KJ/m6MbBo0B/76ZryVNbiD1VS1b2bo4GuY4QCVv/0HqG4fhmniqqpFdflQ12xE8QdBc2KceJfE0behsgbL5cGy8sGNYve1yiUTS/10xS149913p12ez+AGuObOKdmZuXwV8rH27dvH8PAwmUwGj8fD5s2b8Xg83HnnnSiKwqFDh/iLv/gLHnnkEUzTpLbWroje3t7O7bffDsA777zDl770JRwOBw8++CAez8xdvVfrcbbSSQ7OTVIdTrsmSV8nhGqw6uZvWrmg5aG7ib5xAL9fxxn04RlJEnI6iBkG67/wf8z744lbk0rEiV46b+fGWCaYBgoQbL8DX2UYyzBIj1xBr6zFWVlDFpN4fALVqZPLZMilU+B0gWXa+TZOHdZssNs6VFTa51Sxk401R3E3lVi5fD4f3//+9zl06BBf+cpX5v38v/7rv85XvvIVWcpcQV566SXS6fS0Y06nkwcffBBN04jFYiiKQnV1NV1dXSSTSRIJ+4vOqVOnZpxv27ZtfPrTn+bhhx+mpaVl1scsBFClRgKcm6QoClYqTu61V7ASGWr/ave8P8Zd3/3BtMvZ6BD/j8eLw+Od98cStyYRG2N0oIdpXeNVDRSVZGyUVDyGmhgn0duBp7qBzOvPQ7AKpbJh+u6o4iygxw5oqlSs4cvg0FF0N4pDB4du30dRCdbNX0KqWHzbt2/nc5/7HJ2dnTz//PPzvoW7rq5u2kxNJpOhv7+f5uZmKRq4DO3du3fW49lsloMHD1JdXc358+dRVZXvfOc7vPDCC/z5n/85q1bZrWB0XZ+RTOzz+Vi7di2vv/46n/3sZ2f8dy9UPC5F8gq/SZZloVzuxhoaxooOoy7CdK+zKizBzTI1Gumf3ocM7KDFoZMZHSId6ScxGgWXj8Slc5BJQ3ycqta19owNJi6Pj3DrWvSKyuIylBUbhpErMB61Z3UK51VU3IEQvkBoSZ6vmB+JRIK+vj5gcZYJdF2ntbVVgpsVKBKJcObMGQzDIJvNsm/fPmKxGLlcrjgzo+s6mzZtmlaxGOCHP/whzz33HCMjIzPOu3v3/H85Xy5kBucmJaJXsBITOO7cgnr/x8npLvljljGH7iKXSkx2ji90lTfsruFYJoqmguZBaVqD6fZDQysjg31gGliWRTo+jpFJkYvH7CBG1UB3Q3UjSnjKTI2qgQLB6poleKZiPj377LM0Njayc+dOPvvZzy71cMQK8s1vfpO+vj5Wr15NT08PABMTE1y4cIGurq5pt/393/99hoeHqa6unna8qqpqRSUX3yj5TL4Jw52nSA8PMHLsffa9fY4H1t1LVXb2GgOiPISbWhm4cMaeeSkEOYoCmj3bAirkJtfV1bome5kpv6Sl+yoAi8x4YcumPROkBMMQCkM6iaVqdvE3y0T3B6Vx6woWj8d59dVXOXnyJJ2dnQSDwZLuCSTmZrYlpqupqqqiqqoKsFcUXC4Xzc3Nxbo4UzU2Ns5aX+mee+65tQEvczJPeYMGL54jnf8Q6vA0c/jYJU698jbeUNUSj0wsJVVVCdU1oLi99hKTw5nPlXHagYzLDd6KyR+HPlmwz+ECBRRVwx2sorJ1HYrLDZaFomqQTmH1nMOKDoJqBzW6XnoJgeUinU7zyiuvYFkWjz76KOFwmGPHjpV0wTUxN48++ihVVVWEQje+9FxfX8/w8DCtra20tLSwZcuW697ng60bSo3M4NyAVHwCI50E3QMOFzs/8ylCLU207bxrqYcmlgFvoBJvoBKAvvMnZ96gkKNTqGNjmXaQY2TIjOcAi5q2DThdbly+AAPnT4BpgdMNngCKx19M8fFVlGZSYDl4aqiFuQAAEC9JREFU6aWXir8risLXvvY1EonEtJYKojwpisJ9990H2H3J3nvvvTnf99KlS4D9ZWvXrl2AXcSvu7t71tuX8tJUgczg3ACHXkjyVEFVUVWFdRta0ExZnhLTNa7dROPaTdS3bZhSt0Yt9pGyg5x82wWsyf5SlzqYGB1G1TTq120GpxNF11FXrUUJVIKi4K2ql+WpFeyDO1YcDgd9fX3U19cv0YjEclRXV8eePXvYs2fPrD2kZuN0OhkeHuaNN94glUqxZcsW2traZr1tKScXF8gMzg1wOPXJb95YWONRGL0CE6Vb6lrcGlVVqW/fyEDnafuAZaH7KlAdGqZhksukUCyLQHUNo5cvgWkSG+glNthn5++APYsD+erFTnS3a2mejJgXd911F6+++mrx8jPPPMObb77JunXr2LFjWRd4Fkvk3nvv5cCBA0QiEQDC4TAOhwOHw0E8HiedThMOh4lEImSzWcbGxnj55ZdRVfWqy1DlkPMlAc4N8laGSYxE7G/eFTV2cJNflhBiNqqm0bhu83VvZ6ESi/Rh5bKABaZpz+yoKlj5fwGvtOhY0TweD6qqFuvT/OIv/iJNTU189KMfXeKRieVsLgnB2WyWY8eOFXtUmaY5a8Xq+a6YvVxJL6pbkBqPYZo5vEFJMBbzxzAMYpF+kuNjkwctistata3tOFbwEpX0oppkmiYXL16ksbGxJCvJiqUzMDDA+fPnGRsbm/X6PXv2LPKI5s9c30NkBucmGYZBtOM4isstAY6YV5qmUVm/isr6VRjZLIM9XWBkJ693lPbOh3Jy7Ngxent7CQaDEuCIeVVfX099fT2maRKNRtm/f3/xunJYngIJcG5apPMkmDmsZHyphyJKmOZ00njbWgByWTvIURTZG1Aqent7Adi/f/+K/kYtli9VVQmHw8XX18jIyE1tQ1+J5J3yJvmq6uykTyl5LhaJw+lc0UtT4ur8fv9SD0GUicrKSpnBEdcWCNcRCNct9TCEECuYzNoIsXBk+kEIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMARQgghRMmRAEcIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMARQgghRMmRAEcIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMARQgghRMmRAEcIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMARQgghRMmRAEcIcUNM01zqIQghVjDLsrAsa8Efx7HgjyCEKBl9h14Hy8ThD1K7/q6lHo4QYoW5cuUK7777LgB79uxZ0MeSGRwhxJwYuSxgf+sys7mlHYwQYkU6c+bMoj2WzOAIIQDIplOomsbYsf0kDu4D3Q0b76F2yz3obg/ZVBKcHnDo1G+8c6mHK4RYZkzTZGxsDI/Hw+c//3nOnDnDE088wRNPPMHWrVvRdR2n0wnA+vXrF3w8EuAIUcYmRofIJCdwOD1M9F0CLMyT78L4KJgmjF0hcu4oqreCcHMbvvpmPIHgUg9bCLFMpNNpLly4gK7rDAwMMDIywsDAAM899xwAwWCQe+65hxdffJG1a9eydu1aQqEQbW1tCz42CXCEKEPD/d2k+i5CNgO6G0VzgGWBoqKsvRMrVANuL4o/CKoDK5sl0t0JCsRjI8Xz1LS04dRdS/hMhBBLZe/evbMer6+v50/+5E+Ix+Ns3bq1ePz8+fOcP38egM7OTgAqKyv50Ic+tCDjkwBHiDKUHo4ASv6SAqqK0+3FGagiMRZFqagGLFBVULWrnme4v4f61vbFGLIQYpnTNI01a9bQ09PDXXfNbRPCyMgIg4OD1NXVzft4JMARogxVt21kuPMkeEM0btgy7bpQfROWZZGciJHLZpmIXrnqecxslr7OMwRrm/AFAgs9bCHEMuL1ekkkEjz00EP4fL7i8UJ+TTweJxKJcObMGXK5q29MeO+999B1nd27d8/r+JTF2Iu+2Hbs2GEdPHhwqYchRMmIx0YZiwwAoCiKXcPCsijsqgJobN84p3MpinLIsqwdCzHO+SLvIULMr1wux6FDh4hEIle9jcPh4GMf+9h1zzXX9xCZwRFCXJevIoSvIgTA6NAgydgoKAqWaTI1yBFCiNk4HA7uuecewA52XnzxxRlFQ3Vdn9/HnNezCSFKXihcRyg8uV5eirPAQoiF43A4eOyxx6YdM00TVZ3f0nwS4AghbomiKNe/kRBCXMN8BzcglYyFEEIIUYIkwBFCCCFEyZEARwghhBAlRwIcIYQQQpQcCXCEEEIIUXIkwBEi78J//hXOPfUgmcHLSz0UIcQK9Oabb6IoChs3zq3opVhYEuCIsmZZJpf+41OcfvBOTv3vt4me6OHCF59a6mEJIVaQZ555BkVR+PCHPwzA6dOnee2115Z4VEICHFHWYkfewoyO0D2W5LkLUV4fGMexftNSD0sIsYL8xm/8xoxj999//xKMREwlAY4oa+7Vt6PVhAiHA9wedNNe4Sb+2hsY12gMJ4QQUzmdzhnHwuHwEoxETCUBjihrrqpaAp/4d/gDXj7eHGJ1wIWRzDH8P/77Ug9NCLFCJJNJVq1aNe3Y6OjoEo1GFEiAI8pe5UOfYN1zP6fu6SfQK9x4G0PU/Pv/tNTDEkKsEJqm0dPTwwsvvFBsGDmXrthiYUkvKiEAxeGk6b/9JU3/balHIoRYqXbv3k06nV7qYYg8mcERQgghRMmRAEcIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMApEz3/+Ul6/sMnGHjrxaUeihBiBfrUpz6Foiioqko2m13q4QhxXRLglItUCtJpssfe5vKRt5Z6NEKIFebw4cMAWJbFCy+8wIEDB5Z4REJcmwQ4ZSL42/8X7NiFumkHAH2nDtN39tQSj0oIsVJ0d3eza9cu/vAP/xCASCTC3r17SSQSSzwyIWYnAU6ZqNhyL82//VWU6gYU3QWWBbkkfScPLfXQhBArxNtvvz2j79K+ffswDGOJRiTE1UmAU2YaN23HVVm/1MMQQqxQbW1t7NmzZ9oxCXDEciSVjMtQdeMqaFxFJpnA6fYs9XCWlGmaqKrE+ULcqD179pBOpzEMo9ieoFzJ+8jyJAFOGdM93qUewpKxLIu+k4ex+i+AU2fVg59c6iEJseK4XK6lHsKS6u3t5f333wegsbGRbdu2LfGIxFQS4IiSl5qYIDrQA4DD7SGXSoAFGAakJiDnIp1M4CrjgE8IcW2vv/46sVgMr9c7a2J1X1+fBDjLjMypiZIXHewt/p5LFd6YLBSXC6VlI8qqtQz3di3J2IQQy59hGMRiMYBr7hqbmJhYrCGJOZAAR5S0vrNHITkOZm7a8WxHB+PPPguKhqJqkJqg79SRJRqlEGK56u7u5mc/+9mM45Zl8cMf/pCDBw8Wj7366qt0d3cv5vDENUiAI0qbadlb4qeyLOLP/5j4//oBmbOnsderACNLJhlf9CEKIZavVCo16/FoNMo//dM/8b3vfW/a8WPHjmGa5mIMTVyH5OCIklbTvonhy92YufS044HP/Ttc23egb96ClUvDyCC4fQz3XKCmdS0Ol3uJRiyEWE7WrVuHz+fjyJHpM7zV1dX80R/9EXV1dTPuc/z4cbZs2YKiKIs1TDELCXBESXM6ndSvbgNgZLCf5PgIoKCFa/Dc/2EgP3+Ty4DpwcpmiVw4jSdUTaihZcnGLYRYPpqammhqasIwjGnLVXffffest+/p6WFoaIidO3dSUVGxWMMUHyBLVKJsVNY1EKpbNXlAUUBRUHQP1N8GwRp0fwWWaZAYGVq6gQohliVN03j88cdxu68+w6vrOrquk0wm6e/vX8TRiQ+SAEeUFW9glm9TioLicKJoGqrDSahxNeE16xd/cEKIZU9RFDZs2HDV6zOZDHfffTebNm1i7dq1izgy8UES4IiyU9nQjOqYUnnVsoo/qYkYOcMgnUyQmhhfukEKIZatxsZG6urq0DRt1uuPHj2KoiicP3+eTCazyKMTBZKDI8qOxx/A4w9gGDkGO89Mu87Kphnv70FxOMDhILzqNlRNw6GXd8VWIcQkRVHYuXMnAIcOHZqxFDU+Ps6JEycAu4ZOOBwmHA5LO4dFJn9tUbY0zUF9+8bpB8eGYOwKlpGDTIqhC2e4cvEsZi43+0mEEGVt+/bt11yKunDhAu+++y6nTp2S7eOLTGZwRFlLxydAUcEy7aRjfxBMk2m7O02TgY6TODxeqhpacZR5Y0EhxKRcLlescnwtXV1dXLp0iba2Ntavlxy/xSAzOKJsJcZGGRnoBcuwAxzLQnH7UbwV4HCB5gTNQSHayaXTDF++tMSjFkIsF5lMhgMHDjA4ODin21uWRUdHxzXbPYj5IwGOKEujg32MDvQAij2DU6AooKr5f7X876p9O8DIZRiJDJBNpbBkulmIspVIJHjzzTcZGRnB4bixxZBXX32VSCQivasWmCxRibIzePE8RiZfft0y7QCmEOMoyvTWDoXlKzMHWOBwkhyLkhwbAaCmpQ2nLFkJUVaSyST79u0rXs7dYI6eaZocOHAAAFVVeeyxx+Z1fMImMzii7BjZKW0b1HyyTb7oX/H3qYqzOYUtoZPXZ67Sp0YIUbp6e3vn7VySeLxwJMARZcEwckyMjmDmcjg9PjtYmVbDQpn+o6j28hTYAY7mmLw8hcM5ex0MIURpsSyL4eFhBgYGCIfDSz0cMQeyRCVK2vjoEOORyQTA2JBq74KaOklT+N2a5diMXnmTszy6z4/L45vX8QohlpdUKsXhw4eJRqPFY7W1tYDdusEwjFs6/+7du2/p/uLqJMARJWug8wymkV8bn7LslMtepbKowvQg5zpysjwlRMl7+eWXZxwbGbFz8G41uAG7KGB1dfUtn0fMJAGOKFmK5oBCgIOC7gvg8Vfg1HVGBvum5+IU7zTrmT5w0b6sOZ3zOVwhxDLn9Xqpqqpi06ZNdHR0cPHixVvOoZHqxgtHAhxRsupWt2NZFsoHk4aButY2APounoP8G5TT7SGbStq7poqmJB5P2V2lKAqhmvoFG7sQYnnYs2fPrO8jGzZsYMOGDcTjcV555RXADoCSySSWNbepYE3TCIVC8z5mYZMAR5S02YKbqRrXrMOyLCzLRFFUhi5fIpu+ytKTokw7n1P6UwlRFq71PuLz+dizZw+GYaCqKtFolP37988pyHG73dd9jxI3TwIcUfYURUFR7N1QNatWF7+tWaZJLpsFLFBUFAVGIwNkU0l7+UsIIfIKncWrq6t5/PHHiwFOOp0u1snxeDxEIhEOHjwIwLZt25ZmsGVC3qWF+IDCNypFVXG6ps/ShBtblmJIQogVpvA+4na7px2vr69nz549SzGksiPZTUIIIYQoORLgCCGEEKLkSIAjhBBCiJIjAY4QQgghSo4EOEIIIYQoOcpcCxKtJIqiRIBLSz0OIcSsWi3LqlnqQVyLvIcIsazN6T2kJAMcIYQQQpQ3WaISQgghRMmRAEcIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMARQgghRMmRAEcIIYQQJUcCHCGEEEKUHAlwhBBCCFFyJMARQgghRMn5/wFme4UIpWtemAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "score_threshold = 0.25\n",
    "\n",
    "fig, axs = scr.plot_scrublet_results(embedding, \n",
    "                                     scrublet_results['doublet_scores_observed_cells'], \n",
    "                                     scrublet_results['doublet_scores_simulated_doublets'], \n",
    "                                     score_threshold, \n",
    "                                     order_points = True, \n",
    "                                     marker_size = 4)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
